R packages

library(tidyverse)
library(wesanderson)
library(pscl)
library(Matrix)
library(lme4)
library(DHARMa)
library(bbmle) 
library(stats)
library(hrbrthemes)
library(ggpubr)
library(patchwork)
library(FactoMineR)
library(factoextra)
library(ape)
library(vegan)
library(corrplot)
library(ggordiplots)
library(jtools)
library(colorspace)
library(cowplot)
library(ggrepel)
library(ggsci)
library(ggtext)
library(ggthemes)
library(grid)
library(gridExtra)
library(rcartocolor)
library(scico)
library(showtext)
library(sf)
library(spatialreg)
library(spdep)
library(ade4)
library(adespatial)
library(adegraphics)
library(spdep)
library(sp)
library(zetadiv)
library(geosphere)
library(sjPlot)

Raw reads data vizualisation

Relative reads abundance per Class

First, I intended to visualize the number of raw reads obtained for each Class retrieved using BATR01 primers. To do so, I loaded the tabfile I obtained after processing fasta.files with OBITools processing steps and removing the 1st PCR replicate that was not included for sequencing.

batr_raw<- read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/5_Raw_batr_WORep1.txt", header=T, sep="")

batr_raw<-as_tibble(batr_raw)

I needed then to compute the real number of reads obtained in the different PCR replicated without the first one, since sequences were attributed to this replicate due to sequence leaking during sequencing process.

#total read count PCR replicate 1 included
Total_count<-sum(batr_raw$count)
Total_count
## [1] 136704216
#Total read count without PCR replicate 1
Total_count<-sum(batr_raw[,20:1075])
Total_count
## [1] 131883594
#Create a new column with counts corresponding to the total count without rep1
batr_raw<- batr_raw %>% mutate(count_real=rowSums(batr_raw[,20:1075]))

The taxonomic levels at which sequences were assigned using the database, were: Order, family, genus and species. Since the taxonomic level of my interest is the class, I first needed to regroup the data per order and count the number of reads per order.

#determine the number of sequences per order
count_order<- batr_raw %>% group_by(order_name)%>%
  summarise(sum(count_real))

#I exported this output as a file to be able to modify it. 
##write.table(count_order, "C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/pieChart/1_Raw_Batr_count_order.txt", row.names = F)

I attributed Class to each order using Excel by adding a column named “class”. I then loaded this file and computed the total read count for each class.

batr_count_class<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/pieChart/2_Raw_Batr_count_class.txt", header=T, sep="\t")

#count the number of sequences per Class and compute the sequences that it represents: 
df<-batr_count_class %>% group_by(class)%>%
  summarise(sum(sum.count_real.))%>%
  mutate(rra=(`sum(sum.count_real.)`/sum(`sum(sum.count_real.)`))) %>% 
  mutate(rra_fig=c("0.27","0.38","0.02","0.27","0.06","<0.02"))

Barplot of RRA per class

With ggplot2

cbPalette <- c("grey54", "indianred4",rep("grey54",4))

barplot_class <- ggplot(df, aes(x=rra, y=reorder(factor(class, labels=c("Actinopterygii", "Amphibia", "Aves", "Homo sapiens", "Unidentified","Other Mammalia")), rra), fill=class)) +
  geom_bar(stat="identity", width = 0.9) +
  labs(x="RRA") +
  geom_text(aes(label = rra_fig), vjust = 0.45, hjust=-0.30) + #vjust = vertical; hjust = horizontal
  scale_fill_manual(values = cbPalette, name="")+
  theme_ipsum(grid="X") + 
  xlim(c(0,0.5)) + 
  theme(axis.text.x = element_text(size= 11, color= "black", vjust = -1),
        axis.text.y = element_text(size= 11, color= "black"),
        axis.title.x = element_text(size=12, vjust = -3, hjust = 0.5),
        axis.title.y = element_blank(), 
        legend.position = "none")

barplot_class

Relative reads abundance per focal species

Then, I needed to extract the real number of reads for the species of interest and compute their relative reads abundance

sp<-batr_raw %>% filter(best_identity.Batr_3m_final>0.98 & species_name %in% c("Bufo bufo", "Hyla arborea","Lissotriton vulgaris", "Rana temporaria", "Pelophylax ridibundus", "Pelophylax bergeri")) %>% mutate(rra=count_real/sum(count_real))

Since this study was conducted in two different reserves and that we estimated the densities of each of the focal amphibian species per reserve, I was interested in computing the relative read abundance per species and per reserve.

#Yverdon Reserve
yv <- sp %>% select(., contains(c("species_name", "pt1"))) %>% 
  mutate(numreads_yv = rowSums(.[,c(2:342)])) %>%
  mutate(rra = numreads_yv / sum(numreads_yv)) %>% 
  mutate(reserve = "yv")
#Gletterens Reserve
gl <- sp %>% select(., contains(c("species_name","pt2"))) %>% 
  mutate(numreads_gl = rowSums(.[,c(2:331)])) %>%
  mutate(rra = numreads_gl / sum(numreads_gl))%>%
  mutate(reserve="gle")

Barplot of RRA per species and reserve

I plotted the relative read abundance of focal species per reserve. First, I needed to combine the rra from both reserve:

rra_sp<-bind_rows(gl, yv) %>%
  select(-contains("pt"))

#View(rra_sp)

Then, do the barplot per species and reserves:

palette2<-c("#9E6969", "#BD974C", "#BFC66E", "#91B3B8", "#7373B5", "#563B53")

barplot_sp<-ggplot(rra_sp, aes(y=rra, x= factor(reserve, labels = c("Gletterens", "Yverdon")), fill = factor(species_name, levels = c("Bufo bufo", "Hyla arborea", "Lissotriton vulgaris","Pelophylax bergeri", "Pelophylax ridibundus", "Rana temporaria")))) +
  scale_fill_manual(values = palette2, "Species") +
  geom_bar(stat = "identity", color = "black", width = 0.4)+
  ylim(c(0,1)) + ylab("RRA")+ xlab("")+ 
  geom_text(aes(label = round(rra,3)), hjust = -1)+
  theme_ipsum(grid="X") + 
  theme(text= element_text(family ="A"), 
        axis.text.x= element_text(size = 11, color ="black"),
        axis.text.y = element_text(size=11, color="black"),
        axis.title.y = element_text(size=12, hjust= 0.5, vjust = 5),
        axis.title.x = element_blank(),
        legend.position = "none")

barplot_sp

Relation between RRA and estimated densities

Then, we wanted to investigate the relationship between the densities of amphibians and the rra.

Since we cannot differentiate Pelophylax ridibundus from bergeri morphologically speaking and we need to compute the densities of individuals retrieved with the prenuptial migration survey, we pondered the number of Pelophylax sp. individuals by the rra.

#yv ridi: rra= 0.048, berg: rra= 0.111, pelophylax sp. = 1

nridiyv<- 1*(0.048/(0.048+0.111))
nridiyv
## [1] 0.3018868
nbergyv<-1*(0.111/(0.111+0.048))
nbergyv
## [1] 0.6981132
#gle ridi: rra= 0.031, berg.= 0.139, pelophylax sp. = 2

nridigle<-2*(0.031/(0.031+0.139))
nridigle
## [1] 0.3647059
nberggle<-2*(0.139/(0.031+0.139))
nberggle
## [1] 1.635294
density<-read.table("C:/Users/jguenat/Documents/Paper_Master/Data_final/densities_sp_barriers_bergeri_pondere.txt", header=T, dec=",", sep="\t")

To do so the first thing to do is to put rra and densities in a unique file

dens_rra<- left_join(rra_sp, density, by = c("species_name", "reserve"))

Then I performed a linear regression on these term to test the statistical significance of the relationship between the density and the RRA.

mi<-lm(dens_rra$density~dens_rra$rra)
summary(mi)
## 
## Call:
## lm(formula = dens_rra$density ~ dens_rra$rra)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7378 -4.6080 -3.6512  0.1592 29.8268 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)     3.900      3.759   1.038    0.324
## dens_rra$rra   11.638     13.727   0.848    0.416
## 
## Residual standard error: 10.33 on 10 degrees of freedom
## Multiple R-squared:  0.06706,    Adjusted R-squared:  -0.02623 
## F-statistic: 0.7188 on 1 and 10 DF,  p-value: 0.4164
densi_plot<-ggplot(dens_rra, aes(x=density, y=rra))+ 
  geom_point(aes(color=(species_name), shape=factor(reserve, labels = c("Yverdon", "Gletterens"))), size =4, alpha=1)+
  scale_color_manual(values = palette2, "Species")+
  scale_shape_manual(values= c(19,18), "Location")+
  stat_smooth(method="lm", se=F, col="black") + 
  ylim(c(0:1))+ xlim(0,40)+
  theme_minimal()+ ylab("RRA")+ xlab("Density [n/ha]")+
  theme(axis.text = element_text(size=11, color = "black"),
        axis.title.x = element_text(size=12),
        axis.title.y = element_text(size = 12, vjust = -20), 
        axis.ticks = element_blank(),
        legend.text = element_text(size=11, face="italic"),
        legend.title = element_text(size=12),
        legend.position = "none")+
  stat_regline_equation(label.y = 1, label.x =0.01, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 0.90, label.x =0.01, aes(label = ..rr.label..)) + 
  guides(color = guide_legend(order = 1, keyheight = 1.5, keywidth = 1), 
         shape = guide_legend(order = 2, keyheight = 1.5, keywidth = 1))

densi_plot
## `geom_smooth()` using formula = 'y ~ x'

It is also interesting to assess if the RRA is correlated to the number of sampling point where the species were found.

numpoint_plot<-ggplot(dens_rra, aes(x=pres_point, y=rra))+ 
  geom_point(aes(color=(species_name), shape=factor(reserve, labels = c("Yverdon", "Gletterens"))), size =4, alpha=1)+
  scale_color_manual(values = palette2, "Species")+
  scale_shape_manual(values= c(19,18), "Location")+
  stat_smooth(method="lm", se=F, col="black") + 
  ylim(c(0:1))+ xlim(0,25)+
  theme_minimal()+ ylab("RRA")+ xlab("Number of detection points")+
  theme(axis.text = element_text(size=11, color = "black"),
        axis.title.x = element_text(size=12),
        axis.title.y = element_text(size = 12, vjust = 5), 
        axis.ticks = element_blank(),
        legend.text = element_text(size=11, face="italic"),
        legend.title = element_text(size=12),
        legend.position = "right")+
  stat_regline_equation(label.y = 1, label.x =0.01, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 0.90, label.x =0.01, aes(label = ..rr.label..)) + 
  guides(color = guide_legend(order = 1, keyheight = 1.5, keywidth = 1), 
         shape = guide_legend(order = 2, keyheight = 1.5, keywidth = 1))

numpoint_plot
## `geom_smooth()` using formula = 'y ~ x'

Then we statistically determine the relationship between the number of point where the spe. was detected and its rra.

ptdetec<-lm(dens_rra$pres_point~dens_rra$rra)

summary(ptdetec)
## 
## Call:
## lm(formula = dens_rra$pres_point ~ dens_rra$rra)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8919 -3.9481 -0.6081  2.1735 13.0117 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     6.890      1.912   3.604  0.00482 **
## dens_rra$rra    1.661      6.982   0.238  0.81678   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.255 on 10 degrees of freedom
## Multiple R-squared:  0.005627,   Adjusted R-squared:  -0.09381 
## F-statistic: 0.05659 on 1 and 10 DF,  p-value: 0.8168

Figure 2

To produce the fig. 2 of the manuscript, I combined the four plots above.

fig2<-(barplot_class + barplot_sp)/(densi_plot + numpoint_plot) + plot_layout(guides ='collect')# package patchwork

fig2_final<-fig2 + plot_annotation(tag_levels = 'A') & 
  theme(plot.tag = element_text(size = 12), 
        legend.spacing = unit(.5, "cm"), 
        legend.spacing.x = unit(0.5, "cm"))

fig2_final
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Surface covered by the eDNA sampling

Since there is no relationship between RRA and estimated densities of the local amphibian species, I wanted to point out that the sampled surface using eDNA was largely smaller than the total surface of each reserve. Also the eDNA survey was not meant to monitor the density of each species, but rather investigate where these species occur.

So I wanted to compute the surface covered by the eDNA survey: A sampling point represents a circle of 5m diameter. The surface of a point is equal to pi* (rayon)^2

surface <- pi*(2.5)^2

round(surface, 2)
## [1] 19.63

We have 25 sampling point per reserve:

19.63*25# m^2
## [1] 490.75

Environmental data selection

We recorded 13 environmental variable to caracterize the amphibian habitat. However, we only have 50 sampling points and thus for further analyses we need to select the five most relevant variables. To do so, I performed a correlation plot to have an overview of the redunduncy of the variable and then selected among the least correlated variable, the most biologically relevant.

env<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/final_env_Dec20.txt", header=T)
env<-as_tibble(env)

env.corr<-cor(env[,c(5:7, 9:18)])
corrplot(env.corr)

PCA(env[,c(5:7, 9:18)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 13 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

Since there is a lot of variable, I decided to remove the proportion of each vegetation type.

env1.corr<-cor(env[, c(5:7, 9:11, 16:18)])
corrplot(env1.corr)

Then, I’ll check the continuous variables:

env1<-env[, c(5:7, 9:11, 16:18)]
as.tibble(env1)
## # A tibble: 50 × 9
##    X1_submerg X1_emerg X1_emerg_Land mean.water mean.mud forest.dist MeanTemp
##         <int>    <int>         <int>      <dbl>    <dbl>       <dbl>    <dbl>
##  1          4       97            70      2.52     6.17         66.9     16.3
##  2          3       65            15      9.02     4.77         81.6     17.9
##  3          0       95            75      0.856    2.5          32.9     16.7
##  4          0       80            60      2.89     0.838        60.3     17.7
##  5          0       95            55      0.5      1.6         145.      16.8
##  6          1       70            60      2.67     3.37         95.3     19.3
##  7          1       93            35      4.66     1.89         58.1     16.6
##  8          1       60            45      0.856    1.7         140.      20.2
##  9          1       90            10      6.34     4.29         69.0     16.5
## 10          1       75            40      2.74     1.49         60.5     18.6
## # ℹ 40 more rows
## # ℹ 2 more variables: MinTemp <dbl>, MaxTemp <dbl>
continuous<-select_if(env1, is.numeric)
#summary(continuous)# not same scale --> rescaling

#standardize the numeric variables: 

env_scales<-env1 %>% mutate_if(is.numeric, funs(as.numeric(scale(.))))
PCA(env_scales)

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 9 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
par(mfrow=c(3,3))
hist(env_scales$X1_submerg, breaks = 30)
hist(env_scales$X1_emerg, breaks = 30)
hist(env_scales$X1_emerg_Land, breaks = 30)
hist(env_scales$mean.water, breaks = 30)
hist(env_scales$mean.mud, breaks = 30)
hist(env_scales$forest.dist, breaks = 30)
hist(env_scales$MeanTemp, breaks = 30)
hist(env_scales$MinTemp, breaks = 30)
hist(env_scales$MaxTemp, breaks = 30)

par(mfrow=c(3,3))
plot(env_scales$X1_submerg, breaks = 30)
plot(env_scales$X1_emerg, breaks = 30)
plot(env_scales$X1_emerg_Land, breaks = 30)
plot(env_scales$mean.water, breaks = 30)
plot(env_scales$mean.mud, breaks = 30)
plot(env_scales$forest.dist, breaks = 30)
plot(env_scales$MeanTemp, breaks = 30)
plot(env_scales$MinTemp, breaks = 30)
plot(env_scales$MaxTemp, breaks = 30)

pairs(env_scales)

From this I decided to reject Mean water depth, since it is quite correlated to the vegetation as well as the temperature. Moreover it is not necessarily biologically interesting. Indeed, water depth will determine the temperature but also the predation risk (as the more depths the more fishes) but I don’t have any other data to support this, so it is not interesting in my point of view.

I’ll keep the distance to the wintering habitats (estimation of the dispersion ability of sp.) and the mud depth (since important for the wintering habitat of Pelophylax sp. + spots for hunt or hind of sp.).

To determine which other environmental variables I’ll use, I did a last corr plot:

env2.corr<-cor(env_scales)
corrplot(env2.corr)

PCA(env_scales)

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 9 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

I decided to keep: (1) Distance to the wintering habitat; (2) mud depth; (3) mean Water T°C; (4) proportion of emerged vegetation; (5) proportion of submerged vegetation; (6) proportion of emerged land.

env.reduced<-env[, c(1:7, 10:11, 16)]

##write.table(env.reduced,"C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced.txt" ,row.names = F)

env.reduced$points<-as.character(env.reduced$points)
env.scaled<-env.reduced %>% mutate_if(is.numeric, funs(as.numeric(scale(.))))

##write.table(env.scaled,"C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced_scaled.txt" ,row.names = F)

Final Dataset

In order to analyse the effect of environmental factors on the probability of presence of amphibians sp., we first needed to combine the dataset with Pres_Abs of each species with the environmental variables dataset.

env<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced.txt", header=T)
env<-as_tibble(env)

#scaled: 
env2<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced_scaled.txt", header=T)
env2<-as_tibble(env2)

batr<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/8_Batr_presabs.txt", header=T,sep="\t", check.names = F)
batr<-as_tibble(batr)

batr1<-batr %>% select(species_name, 11:60) %>%
  pivot_longer(2:51, names_to="points", values_to="presabs")

env$points<-as.character(env$points)
evabund<-left_join(env, batr1, by="points")

#scaled: 
env2$points<-as.character(env2$points)
evabund2<-left_join(env2, batr1, by="points")

#For convenience for further analyses, I will put each species as a column instead of in line: 

dat1<-evabund2

bbufo<- dat1[dat1$species_name=="Bufo bufo",]
b<-bbufo %>% select(points, presabs) %>%
  rename(Bbuf_pres=presabs)

dat1<-right_join(dat1, b)

Harbo<- dat1[dat1$species_name=="Hyla arborea",]
h<-Harbo %>% select(points, presabs) %>%
  rename(Harbo_pres=presabs)

dat1<-right_join(dat1, h)

Lvulga<- dat1[dat1$species_name=="Lissotriton vulgaris",]
L<-Lvulga %>% select(points, presabs) %>% 
  rename(Lvulga_pres=presabs)

dat1<-right_join(dat1, L)

Pridi<- dat1[dat1$species_name=="Pelophylax ridibundus",]
P<-Pridi %>% select(points, presabs) %>% 
  rename(Pridi_pres=presabs)
dat1<-right_join(dat1, P)

Pberg<- dat1[dat1$species_name=="Pelophylax bergeri",]
pb<-Pberg %>% select(points, presabs) %>% 
  rename(Pberg_pres=presabs)
dat1<-right_join(dat1, pb)

Rtempo<- dat1[dat1$species_name=="Rana temporaria",]
R<- Rtempo %>% select(points, presabs) %>% 
  rename(Rtemp_pres=presabs)
dat1<-right_join(dat1, R)
#View(dat1)

x<- dat1 %>% select(!c(species_name, presabs)) %>% distinct()
#View(x)

##write.table(x,"C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/2022_FinalEnv_Batr_scaled.txt", row.names = F)

Statistical analyses

data<- read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/2022_FinalEnv_Batr_scaled.txt", header=T, sep="")

Zeta diversity

“Zeta diversity is a concept in ecology used to understand the turnover of species across multiple sites or sampling units. Unlike traditional diversity measures that typically focus on individual sites (alpha diversity) or pairs of sites (beta diversity), zeta diversity generalizes the concept to any number of sites, providing a more holistic view of species turnover across landscapes or communities.

Definition: Zeta diversity (ζ-diversity) measures the number of shared species among a given number of sites. For example, ζ₁ represents the total number of species observed in any single site, ζ₂ the number of species shared by any two sites, ζ₃ the number shared by any three sites, and so on.

Purpose: It provides insights into the scaling of biodiversity, showing how species turnover changes with the number of sites considered. It can help in understanding the spatial distribution of biodiversity, the effects of habitat fragmentation, and species co-occurrence patterns.

Calculation:

ζ₁ is simply the average species richness across sites. ζ₂ is the average number of species shared between any pair of sites. ζ₃ is the average number of species shared among any three sites, and so on. The values can be calculated for k sites (ζₖ) and usually decrease as k increases, reflecting the likelihood that more sites share fewer species.

Zeta diversity is implemented in R using the zetadiv package.

Then, we need to prepare the data to compute the zeta diversity. For that we need one matrix with the species pres/abs and a matrix with the geographical distances.”

Yverdon

select the data:

data2<- read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/2024_coordadded.txt", header=T, sep="")


data_1 <- data2 %>% filter(Location == "Cheseaux")
xy = data_1[, c("x", "y")]
xy_coord <- st_as_sf(xy, coords = c("x", "y"), crs = 4326)
coordinates_utm <- st_transform(xy_coord, crs = 32632)
xy <- as.data.frame(st_coordinates(coordinates_utm))
occurences = data_1[, 11:16]

for each order (2, 3, 5, 10), a linear regression and a gam are computed by default distance = Euclidean.

set.seed(1)
zeta.ddecay.lm.fine2 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 2, confint.level = 0.95, method = "mean", plot = T)

set.seed(1)
zeta.ddecay.lm.fine3 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 3, confint.level = 0.95, method = "mean", plot = T)

set.seed(1)
zeta.ddecay.lm.fine5 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 5, confint.level = 0.95, method = "mean", plot = T)

set.seed(1)
zeta.ddecay.lm.fine10 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 10, confint.level = 0.95, method = "mean", plot = T)

Zeta 2

To be able to modify the graphs, we need to extract the information:

distance_z2 <- zeta.ddecay.lm.fine2$distance
zeta2<- zeta.ddecay.lm.fine2$zeta.val

data_plot_z2_yv <- data.frame(distance = distance_z2, zeta2= zeta2)

coeff_z2_yv <- coef(zeta.ddecay.lm.fine2$reg)
intercept_z2_yv <- coeff_z2_yv[1]
slope_z2_yv <- coeff_z2_yv[2]

#Extract predicted values and standar errors: 
predictions_z2_yv <- predict(zeta.ddecay.lm.fine2$reg, newdata = data_plot_z2_yv, se.fit = TRUE)

# Extract predicted values and standard errors
data_plot_z2_yv$predicted <- predictions_z2_yv$fit
data_plot_z2_yv$se <- predictions_z2_yv$se.fit

# Calculate confidence intervals
data_plot_z2_yv$lower <- data_plot_z2_yv$predicted - 1.96 * data_plot_z2_yv$se
data_plot_z2_yv$upper <- data_plot_z2_yv$predicted + 1.96 * data_plot_z2_yv$se

Graphical representation:

zeta2_yv <- ggplot(data_plot_z2_yv, aes(x = distance, y = zeta2)) +
  geom_point(color = "black", size = 2, alpha = 0.25) + 
  geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") + 
  labs(x = "Distance between two sites (m)", 
       y = expression("Number of species shared by two sites (ζ"[2]*")"), 
       title = "Yverdon") +
  theme_light() +                                 
  theme(axis.title.x = element_text(size = 11), 
    axis.title.y = element_text(size = 11),         
    plot.title = element_text(hjust = 0.5, size = 12))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Zeta 3

distance_z3 <- zeta.ddecay.lm.fine3$distance
zeta3<- zeta.ddecay.lm.fine3$zeta.val

data_plot_z3_yv <- data.frame(distance = distance_z3, zeta3= zeta3)

coeff_z3_yv <- coef(zeta.ddecay.lm.fine3$reg)
intercept_z3_yv <- coeff_z3_yv[1]
slope_z3_yv <- coeff_z3_yv[2]

#Extract predicted values and standar errors: 
predictions_z3_yv <- predict(zeta.ddecay.lm.fine3$reg, newdata = data_plot_z3_yv, se.fit = TRUE)

# Extract predicted values and standard errors
data_plot_z3_yv$predicted <- predictions_z3_yv$fit
data_plot_z3_yv$se <- predictions_z3_yv$se.fit

# Calculate confidence intervals
data_plot_z3_yv$lower <- data_plot_z3_yv$predicted - 1.96 * data_plot_z3_yv$se
data_plot_z3_yv$upper <- data_plot_z3_yv$predicted + 1.96 * data_plot_z3_yv$se

zeta3_yv <- ggplot(data_plot_z3_yv, aes(x = distance, y = zeta3)) +
  geom_point(color = "black", size = 2, alpha = 0.25) + 
  geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") + 
  labs(x = "Mean distance between three sites (m)", 
       y = expression("Number of species shared by three sites (ζ"[3]*")"), 
       title = "Yverdon") +
  theme_light() +                                 
  theme(axis.title.x = element_text(size = 11), 
        axis.title.y = element_text(size = 11),         
        plot.title = element_text(hjust = 0.5, size = 12))

Zeta 5

distance_z5 <- zeta.ddecay.lm.fine5$distance
zeta5<- zeta.ddecay.lm.fine5$zeta.val

data_plot_z5_yv <- data.frame(distance = distance_z5, zeta5= zeta5)

coeff_z5_yv <- coef(zeta.ddecay.lm.fine5$reg)
intercept_z5_yv <- coeff_z5_yv[1]
slope_z5_yv <- coeff_z5_yv[2]

#Extract predicted values and standar errors: 
predictions_z5_yv <- predict(zeta.ddecay.lm.fine5$reg, newdata = data_plot_z5_yv, se.fit = TRUE)

# Extract predicted values and standard errors
data_plot_z5_yv$predicted <- predictions_z5_yv$fit
data_plot_z5_yv$se <- predictions_z5_yv$se.fit

# Calculate confidence intervals
data_plot_z5_yv$lower <- data_plot_z5_yv$predicted - 1.96 * data_plot_z5_yv$se
data_plot_z5_yv$upper <- data_plot_z5_yv$predicted + 1.96 * data_plot_z5_yv$se

zeta5_yv <- ggplot(data_plot_z5_yv, aes(x = distance, y = zeta5)) +
  geom_point(color = "black", size = 2, alpha = 0.25) + 
  geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") + 
  labs(x = "Mean distance between five sites (m)", 
       y = expression("Number of species shared by five sites (ζ"[5]*")"), 
       title = "Yverdon") +
  theme_light() +                                 
  theme(axis.title.x = element_text(size = 11), 
        axis.title.y = element_text(size = 11),         
        plot.title = element_text(hjust = 0.5, size = 12))

Gletterens

data_2 <- data2 %>% filter(Location == "Glettrens")
xy = data_2[, c("x", "y")]
xy_coord <- st_as_sf(xy, coords = c("x", "y"), crs = 4326)
coordinates_utm <- st_transform(xy_coord, crs = 32632)
xy <- as.data.frame(st_coordinates(coordinates_utm))
occurences = data_2[, 11:16]

for each order (2, 3, 5, 10), a linear regression and a gam are computed by default distance = Euclidean

set.seed(1)
zeta.ddecay.lm.fine2 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 2, confint.level = 0.95, method = "mean", plot = T)

set.seed(1)
zeta.ddecay.lm.fine3 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 3, confint.level = 0.95, method = "mean", plot = T)

set.seed(1)
zeta.ddecay.lm.fine5 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 5, confint.level = 0.95, method = "mean", plot = T)

set.seed(1)
zeta.ddecay.lm.fine10 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 10, confint.level = 0.95, method = "mean", plot = T)

Zeta 2

distance_z2 <- zeta.ddecay.lm.fine2$distance
zeta2<- zeta.ddecay.lm.fine2$zeta.val

data_plot_z2_gle <- data.frame(distance = distance_z2, zeta2= zeta2)

coeff_z2_gle <- coef(zeta.ddecay.lm.fine2$reg)
intercept_z2_gle <- coeff_z2_gle[1]
slope_z2_gle <- coeff_z2_gle[2]

#Extract predicted values and standar errors: 
predictions_z2_gle <- predict(zeta.ddecay.lm.fine2$reg, newdata = data_plot_z2_gle, se.fit = TRUE)

# Extract predicted values and standard errors
data_plot_z2_gle$predicted <- predictions_z2_gle$fit
data_plot_z2_gle$se <- predictions_z2_gle$se.fit

# Calculate confidence intervals
data_plot_z2_gle$lower <- data_plot_z2_gle$predicted - 1.96 * data_plot_z2_gle$se
data_plot_z2_gle$upper <- data_plot_z2_gle$predicted + 1.96 * data_plot_z2_gle$se

zeta2_gle <- ggplot(data_plot_z2_gle, aes(x = distance, y = zeta2)) +
  geom_point(color = "black", size = 2, alpha = 0.25) + 
  geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") + 
  labs(x = "Distance between two sites (m)", 
       y = expression("Number of species shared by two sites (ζ"[2]*")"), 
       title = "Gletterens") +
  theme_light() +                                 
  theme(axis.title.x = element_text(size = 11), 
        axis.title.y = element_text(size = 11),         
        plot.title = element_text(hjust = 0.5, size = 12))

Zeta 3

distance_z3 <- zeta.ddecay.lm.fine3$distance
zeta3<- zeta.ddecay.lm.fine3$zeta.val

data_plot_z3_gle <- data.frame(distance = distance_z3, zeta3= zeta3)

coeff_z3_gle <- coef(zeta.ddecay.lm.fine3$reg)
intercept_z3_gle <- coeff_z3_gle[1]
slope_z3_gle <- coeff_z3_gle[2]

#Extract predicted values and standar errors: 
predictions_z3_gle <- predict(zeta.ddecay.lm.fine3$reg, newdata = data_plot_z3_gle, se.fit = TRUE)

# Extract predicted values and standard errors
data_plot_z3_gle$predicted <- predictions_z3_gle$fit
data_plot_z3_gle$se <- predictions_z3_gle$se.fit

# Calculate confidence intervals
data_plot_z3_gle$lower <- data_plot_z3_gle$predicted - 1.96 * data_plot_z3_gle$se
data_plot_z3_gle$upper <- data_plot_z3_gle$predicted + 1.96 * data_plot_z3_gle$se

zeta3_gle <- ggplot(data_plot_z3_gle, aes(x = distance, y = zeta3)) +
  geom_point(color = "black", size = 2, alpha = 0.25) + 
  geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") + 
  labs(x = "Mean distance between three sites (m)", 
       y = expression("Number of species shared by three sites (ζ"[3]*")"), 
       title = "Gleterrens") +
  theme_light() +                                 
  theme(axis.title.x = element_text(size = 11), 
        axis.title.y = element_text(size = 11),         
        plot.title = element_text(hjust = 0.5, size = 12))

Zeta 5

distance_z5 <- zeta.ddecay.lm.fine5$distance
zeta5<- zeta.ddecay.lm.fine5$zeta.val

data_plot_z5_gle <- data.frame(distance = distance_z5, zeta5= zeta5)

coeff_z5_gle <- coef(zeta.ddecay.lm.fine5$reg)
intercept_z5_gle <- coeff_z5_gle[1]
slope_z5_gle <- coeff_z5_gle[2]

#Extract predicted values and standar errors: 
predictions_z5_gle <- predict(zeta.ddecay.lm.fine5$reg, newdata = data_plot_z5_gle, se.fit = TRUE)

# Extract predicted values and standard errors
data_plot_z5_gle$predicted <- predictions_z5_gle$fit
data_plot_z5_gle$se <- predictions_z5_gle$se.fit

# Calculate confidence intervals
data_plot_z5_gle$lower <- data_plot_z5_gle$predicted - 1.96 * data_plot_z5_gle$se
data_plot_z5_gle$upper <- data_plot_z5_gle$predicted + 1.96 * data_plot_z5_gle$se

zeta5_gle <- ggplot(data_plot_z5_gle, aes(x = distance, y = zeta5)) +
  geom_point(color = "black", size = 2, alpha = 0.25) + 
  geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") + 
  labs(x = "Mean distance between five sites (m)", 
       y = expression("Number of species shared by five sites (ζ"[5]*")"), 
       title = "Gletterens") +
  theme_light() +                                 
  theme(axis.title.x = element_text(size = 11), 
        axis.title.y = element_text(size = 11),         
        plot.title = element_text(hjust = 0.5, size = 12))

Plot zeta diversity

Fig_S4<-(zeta2_yv | zeta3_yv | zeta5_yv) / (zeta2_gle | zeta3_gle | zeta5_gle)

Fig_S4 + plot_annotation(tag_levels = 'A')

db-RDA replacing 0 by 0.001

We asked if the species assembly could be explained by the habitat characteristics proxied by the environmental variables we measured. We performed a distance-based redundancy analysis.

I wanted to have a first look at how the proportion of emerged, submerged vegetation and emerged land affect the species assembly, to determine if these variables are all important ones.

dat<-data

#dat[dat == 0]<- 0.001

sp=dat[,11:16]
sp001<- sp + 0.001
env=dat[,5:10]

#sp.dist001<-vegdist(sp001, method = "jaccard")

colnames(sp)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(sp001)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(env)<-c("SV", "EV", "EL", "MD", "FD", "WT")

rankindex(env, sp001, indices = c("jaccard", "bray"), stepacross= FALSE, method = "spearman")
##     jaccard        bray 
## -0.03328199 -0.03293306
#Capscale function from vegan package 

dbRDA = capscale(sp001 ~ SV + EV + EL + MD + FD + WT,
                 env, dist="jaccard")

summary(dbRDA)
## 
## Call:
## capscale(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env,      distance = "jaccard") 
## 
## Partitioning of squared Jaccard distance:
##               Inertia Proportion
## Total          16.671      1.000
## Constrained     2.268      0.136
## Unconstrained  14.403      0.864
## 
## Eigenvalues, and their contribution to the squared Jaccard distance 
## 
## Importance of components:
##                          CAP1    CAP2    CAP3    CAP4     CAP5      CAP6   MDS1
## Eigenvalue            0.99265 0.64674 0.35951 0.24218 0.023288 0.0031386 5.7253
## Proportion Explained  0.05955 0.03880 0.02157 0.01453 0.001397 0.0001883 0.3434
## Cumulative Proportion 0.05955 0.09834 0.11991 0.13443 0.135831 0.1360188 0.4795
##                         MDS2   MDS3   MDS4    MDS5    MDS6   MDS7      MDS8
## Eigenvalue            3.0656 2.2314 1.8295 0.71615 0.54941 0.2701 0.0155887
## Proportion Explained  0.1839 0.1339 0.1097 0.04296 0.03296 0.0162 0.0009351
## Cumulative Proportion 0.6634 0.7972 0.9069 0.94991 0.98286 0.9991 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CAP1   CAP2   CAP3   CAP4    CAP5     CAP6
## Eigenvalue            0.9926 0.6467 0.3595 0.2422 0.02329 0.003139
## Proportion Explained  0.4378 0.2852 0.1586 0.1068 0.01027 0.001384
## Cumulative Proportion 0.4378 0.7230 0.8815 0.9883 0.99862 1.000000
RsquareAdj(dbRDA)
## $r.squared
## [1] 0.1360188
## 
## $adj.r.squared
## [1] 0.01546327

Since the relation between the probability of presences of species and the environmental factors composing their habitat might not be linear (i.e. species might have an unimodal distribution, with an optimum and some tolerance to the variation), we tested the effect on the model of the addition of quadratic terms of all environmental variables. We used the polyvar() function proposed in “Numerical Ecology with R” from Bocard, Gillet & Legendre (2018).

# adding the quadratic effect of environmental variables: 
env.square<- polyvars(env, degr = 2)

rankindex(env.square, sp001, indices = c("jaccard", "bray"), stepacross= FALSE, method = "spearman")
##     jaccard        bray 
## -0.04761048 -0.04725647
dbRDA.square<- capscale(sp001 ~ ., env.square, distance = "jaccard")
summary(dbRDA.square)
## 
## Call:
## capscale(formula = sp001 ~ SV.1 + SV.2 + EV.1 + EV.2 + EL.1 +      EL.2 + MD.1 + MD.2 + FD.1 + FD.2 + WT.1 + WT.2, data = env.square,      distance = "jaccard") 
## 
## Partitioning of squared Jaccard distance:
##               Inertia Proportion
## Total          16.671      1.000
## Constrained     4.568      0.274
## Unconstrained  12.102      0.726
## 
## Eigenvalues, and their contribution to the squared Jaccard distance 
## 
## Importance of components:
##                         CAP1    CAP2    CAP3    CAP4     CAP5     CAP6
## Eigenvalue            2.1291 1.33975 0.52667 0.39283 0.123676 0.038943
## Proportion Explained  0.1277 0.08037 0.03159 0.02356 0.007419 0.002336
## Cumulative Proportion 0.1277 0.20808 0.23967 0.26324 0.270657 0.272993
##                            CAP7      CAP8   MDS1   MDS2   MDS3    MDS4    MDS5
## Eigenvalue            0.0155790 0.0018441 4.5776 2.8056 1.8136 1.45259 0.65768
## Proportion Explained  0.0009345 0.0001106 0.2746 0.1683 0.1088 0.08714 0.03945
## Cumulative Proportion 0.2739273 0.2740379 0.5486 0.7169 0.8257 0.91286 0.95231
##                          MDS6   MDS7      MDS8
## Eigenvalue            0.53390 0.2468 0.0143806
## Proportion Explained  0.03203 0.0148 0.0008626
## Cumulative Proportion 0.98433 0.9991 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                        CAP1   CAP2   CAP3    CAP4    CAP5     CAP6    CAP7
## Eigenvalue            2.129 1.3398 0.5267 0.39283 0.12368 0.038943 0.01558
## Proportion Explained  0.466 0.2933 0.1153 0.08599 0.02707 0.008524 0.00341
## Cumulative Proportion 0.466 0.7593 0.8746 0.96059 0.98766 0.996186 0.99960
##                            CAP8
## Eigenvalue            0.0018441
## Proportion Explained  0.0004037
## Cumulative Proportion 1.0000000
RsquareAdj(dbRDA.square)
## $r.squared
## [1] 0.2740379
## 
## $adj.r.squared
## [1] 0.03859079

The adjusted R.squared is way better with the quadratic effect of the environmental variables, thus we decided to keep them.

Ordination plot of db-RDA

First, the triplot for the db-RDA without quadratic effect:

shapes = c(22, 23, 24) 
shapes <- shapes[as.factor(dat$environment)]

colors <- c("red", "darkgreen", "blue")
colors <- colors[as.factor(dat$environment)]

Names = paste("N", 1:4, sep=" = ")

plot(dbRDA, type= "n", scaling = 3, las=1)
points(dbRDA, display = "species", pch=20, cex = 1.5, col="#ff7f00", scaling=3)
text(dbRDA, display = "species", scaling= 3, labels = colnames(sp), cex=0.8, pos= 3)
points(dbRDA, display = "sites", pch= shapes, cex= 1.2, col= colors, scaling=3)
text(dbRDA, scaling= 3, display="bp", col="#0868ac", cex= 0.8)

"#0868ac"
## [1] "#0868ac"
# Legend
legend("topright", legend = levels(as.factor(dat$environment)),
 col = c("red", "darkgreen", "blue"),
 pch = c(22, 23, 24) )

Significance of db-RDA

We tested whether the influences of the environmental variables on the species assemblage were significant by performing a permanova (i.e., permutational multivariate analysis of variance) using the function adonis2 from the package ‘vegan’.

#PERMANOVA!!!! 
adonis2(sp001 ~ SV + EV + EL + MD + FD + WT,
        env, dist="jaccard", permutations = 999)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env, permutations = 999, dist = "jaccard")
##          Df SumOfSqs      R2      F Pr(>F)
## Model     6   1.9105 0.13564 1.1246  0.355
## Residual 43  12.1746 0.86436              
## Total    49  14.0851 1.00000
anova(dbRDA, permutations=how(nperm = 999))
## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env, distance = "jaccard")
##          Df SumOfSqs      F Pr(>F)
## Model     6   2.2675 1.1283  0.295
## Residual 43  14.4030
#with quadratic effect
anova(dbRDA.square, permutations=how(nperm = 999))
## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = sp001 ~ SV.1 + SV.2 + EV.1 + EV.2 + EL.1 + EL.2 + MD.1 + MD.2 + FD.1 + FD.2 + WT.1 + WT.2, data = env.square, distance = "jaccard")
##          Df SumOfSqs      F Pr(>F)
## Model    12   4.5684 1.1639  0.226
## Residual 37  12.1022
adonis2(sp001 ~ .,
        env.square, dist="jaccard", permutations = 999)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = sp001 ~ ., data = env.square, permutations = 999, dist = "jaccard")
##          Df SumOfSqs      R2      F Pr(>F)
## Model    12   3.9623 0.28131 1.2069  0.236
## Residual 37  10.1228 0.71869              
## Total    49  14.0851 1.00000

##RDA

dat<-data

sp=dat[,11:16]
sp001<- sp + 0.001
env=dat[,5:10]

colnames(sp)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(sp001)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(env)<-c("SV", "EV", "EL", "MD", "FD", "WT")

#RDA function from vegan package 

RDA = rda(sp001 ~ SV + EV + EL + MD + FD + WT,
                 env)
summary(RDA)
## 
## Call:
## rda(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          1.0898     1.0000
## Constrained    0.1471     0.1349
## Unconstrained  0.9427     0.8651
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          RDA1    RDA2    RDA3    RDA4      RDA5      RDA6
## Eigenvalue            0.06123 0.03067 0.02875 0.02468 0.0009816 0.0007378
## Proportion Explained  0.05619 0.02814 0.02639 0.02264 0.0009008 0.0006770
## Cumulative Proportion 0.05619 0.08433 0.11072 0.13336 0.1342601 0.1349371
##                          PC1    PC2    PC3     PC4     PC5     PC6
## Eigenvalue            0.3472 0.2523 0.1616 0.08730 0.05446 0.03985
## Proportion Explained  0.3186 0.2315 0.1483 0.08011 0.04997 0.03657
## Cumulative Proportion 0.4536 0.6850 0.8333 0.91346 0.96343 1.00000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                          RDA1    RDA2    RDA3    RDA4      RDA5      RDA6
## Eigenvalue            0.06123 0.03067 0.02875 0.02468 0.0009816 0.0007378
## Proportion Explained  0.41640 0.20856 0.19554 0.16781 0.0066754 0.0050175
## Cumulative Proportion 0.41640 0.62496 0.82050 0.98831 0.9949825 1.0000000
RsquareAdj(RDA)
## $r.squared
## [1] 0.1349371
## 
## $adj.r.squared
## [1] 0.01423067
# adding the quadratic effect of environmental variables: 

env.square<- polyvars(env, degr = 2)
RDA.square<- rda(sp001 ~ ., env.square)
summary(RDA.square)
## 
## Call:
## rda(formula = sp001 ~ SV.1 + SV.2 + EV.1 + EV.2 + EL.1 + EL.2 +      MD.1 + MD.2 + FD.1 + FD.2 + WT.1 + WT.2, data = env.square) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          1.0898      1.000
## Constrained    0.2648      0.243
## Unconstrained  0.8250      0.757
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                         RDA1    RDA2    RDA3    RDA4    RDA5     RDA6    PC1
## Eigenvalue            0.1119 0.06855 0.03992 0.02988 0.01117 0.003431 0.2891
## Proportion Explained  0.1027 0.06290 0.03663 0.02741 0.01025 0.003148 0.2653
## Cumulative Proportion 0.1027 0.16558 0.20221 0.22962 0.23987 0.243021 0.5083
##                          PC2    PC3     PC4     PC5     PC6
## Eigenvalue            0.2232 0.1507 0.07658 0.05212 0.03323
## Proportion Explained  0.2048 0.1383 0.07027 0.04782 0.03049
## Cumulative Proportion 0.7132 0.8514 0.92169 0.96951 1.00000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1    RDA2    RDA3    RDA4    RDA5     RDA6
## Eigenvalue            0.1119 0.06855 0.03992 0.02988 0.01117 0.003431
## Proportion Explained  0.4225 0.25884 0.15071 0.11280 0.04218 0.012956
## Cumulative Proportion 0.4225 0.68135 0.83206 0.94487 0.98704 1.000000
RsquareAdj(RDA.square)
## $r.squared
## [1] 0.2430214
## 
## $adj.r.squared
## [1] -0.002485127
# with squared effect: 

Names = paste("N", 1:4, sep=" = ")

plot(dbRDA.square, type= "n", scaling = 2, las=1, xlab = "CAP 1 (12.77%)", ylab= "CAP 2 (8.04%)", ylim = c(-2.5, 2.5), xlim = c(-4, 4))
points(dbRDA.square, display = "species", pch=19, cex = 1.5, col=palette2, scaling=3)
text(dbRDA.square, display = "species", scaling= 3, labels = colnames(sp), cex=0.8, pos= 2)
points(dbRDA.square, display = "sites", pch = 20, cex= 1, col= "black", scaling=3)
text(dbRDA.square, scaling= 3, display="bp", col="#0868ac", cex= 0.8)

Generalized Linear Models

We aimed at characterizing the breeding habitat of each of the local amphibian species. To this end, we tested the effect of the recorded environmental variables on the probability of presence of each amphibian species. We performed separate binomial generalized linear models (GLMs) for each species, using the presence-absence of the focal species as response variable and the environmental variables as explanatory variables.

based on the graphs chunk 45, we can observe that the distribution of points does not have a linear pattern. Thus, I’ll add the quadratic effect of these variables in the following models.

Bufo bufo

length(data$points[data$Bbuf_pres==1])
## [1] 19

We wanted to analyze the effect of other amphibians species on the probability of presence of the focal sp. Since there is 5 other species and that we are already testing 6 environmental variable, we would end with 11 explanatory variable that is way too much for 50 sampling points. What we did then, is to first do a PCA on the presence of other amphibian sp and then use this as an explanatory variable in the model.

pcaB<-prcomp(data[,c(12:16)])
pcaB
## Standard deviations (1, .., p=5):
## [1] 0.6281888 0.4309699 0.3815493 0.2621920 0.2338955
## 
## Rotation (n x k) = (5 x 5):
##                    PC1         PC2          PC3        PC4         PC5
## Harbo_pres   0.1419059  0.10666751 -0.200815801  0.7733082 -0.57458873
## Lvulga_pres -0.1818155 -0.95494934 -0.006808097  0.2211551  0.07783914
## Pridi_pres  -0.1483413  0.15594608 -0.672260323  0.3058557  0.63890062
## Pberg_pres  -0.6728335  0.01282751 -0.453510094 -0.3459734 -0.47091568
## Rtemp_pres   0.6870923 -0.22849523 -0.549564661 -0.3739513 -0.18393825
str(pcaB)
## List of 5
##  $ sdev    : num [1:5] 0.628 0.431 0.382 0.262 0.234
##  $ rotation: num [1:5, 1:5] 0.142 -0.182 -0.148 -0.673 0.687 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "Harbo_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:5] 0.08 0.24 0.12 0.34 0.56
##   ..- attr(*, "names")= chr [1:5] "Harbo_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:50, 1:5] -0.927 -0.779 0.581 0.723 -0.106 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
PCA(data[,c(12:16)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
pcaB$x[,1]#is what we use as explanatory variable. 
##  [1] -0.9270989 -0.7787576  0.5811682  0.7230741 -0.1059241 -0.1059241
##  [7] -0.1059241  0.5811682 -0.9605731 -0.2877396 -0.7787576  0.7230741
## [13]  0.5811682 -0.9605731 -0.9270989 -0.1059241 -0.7787576 -0.1059241
## [19]  0.3993527 -0.7787576 -0.1059241  0.5811682 -0.9605731 -0.9270989
## [25]  0.5811682  0.5811682  0.3993527  0.3993527 -0.1059241  0.5811682
## [31] -0.0916653 -0.7787576  0.5811682  0.5747328 -0.4218221  0.5811682
## [37]  0.5811682  0.5811682 -0.9605731  0.5811682 -0.7787576 -0.2400066
## [43]  0.5811682  0.7230741 -0.9605731  0.5811682  0.3993527  0.5811682
## [49]  0.5811682  0.3993527
#with PC1 as explanatory variable: 

Bbsp<-glm(data$Bbuf_pres ~ pcaB$x[,1], family = "binomial")
testZeroInflation(simulateResiduals(Bbsp))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99871, p-value = 1
## alternative hypothesis: two.sided
simulateResiduals(Bbsp, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2562227 0.2492455 0.7444733 0.9873053 0.03247239 0.0819542 0.4597335 0.06559149 0.0448016 0.1499663 0.4141015 0.2128683 0.5981021 0.7153422 0.9464596 0.4942223 0.5403071 0.1554513 0.9880159 0.2774713 ...
summary(Bbsp)
## 
## Call:
## glm(formula = data$Bbuf_pres ~ pcaB$x[, 1], family = "binomial")
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.4910     0.2919  -1.682   0.0925 .
## pcaB$x[, 1]   0.1775     0.4723   0.376   0.7070  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 66.406  on 49  degrees of freedom
## Residual deviance: 66.264  on 48  degrees of freedom
## AIC: 70.264
## 
## Number of Fisher Scoring iterations: 4
Bb<-glm(data$Bbuf_pres~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land + data$mean.mud + data$forest.dist + data$MeanTemp + pcaB$x[,1], family = "binomial")
testZeroInflation(simulateResiduals(Bb))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99512, p-value = 1
## alternative hypothesis: two.sided
simulateResiduals(Bb, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2082843 0.2280974 0.7469069 0.9897005 0.03628088 0.06526891 0.5528821 0.06559149 0.04862612 0.1682549 0.426804 0.1978776 0.5380052 0.7645423 0.9301942 0.5498223 0.5339506 0.09838691 0.9897279 0.2477422 ...
drop1(Bb)# PCO1 --> other sp. does not explain the probability of presence of B. bufo. 
## Single term deletions
## 
## Model:
## data$Bbuf_pres ~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land + 
##     data$mean.mud + data$forest.dist + data$MeanTemp + pcaB$x[, 
##     1]
##                    Df Deviance    AIC
## <none>                  60.484 76.484
## data$X1_submerg     1   60.762 74.762
## data$X1_emerg       1   61.076 75.076
## data$X1_emerg_Land  1   61.347 75.347
## data$mean.mud       1   63.448 77.448
## data$forest.dist    1   60.637 74.637
## data$MeanTemp       1   60.518 74.518
## pcaB$x[, 1]         1   60.491 74.491
Bb<-glm(Bbuf_pres~ X1_submerg + X1_emerg + X1_emerg_Land +
          mean.mud + forest.dist + MeanTemp, data, family = "binomial")
testZeroInflation(simulateResiduals(Bb))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99538, p-value = 1
## alternative hypothesis: two.sided
simulateResiduals(Bb, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2082843 0.2235656 0.7493404 0.9898203 0.03628088 0.06526891 0.5588917 0.06643783 0.04835294 0.1682549 0.421723 0.1978776 0.5530294 0.761028 0.9261278 0.5498223 0.5307723 0.09346756 0.9898502 0.2427874 ...
drop1(Bb)# Mean_Temp
## Single term deletions
## 
## Model:
## Bbuf_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + 
##     forest.dist + MeanTemp
##               Df Deviance    AIC
## <none>             60.491 74.491
## X1_submerg     1   60.777 72.777
## X1_emerg       1   61.134 73.134
## X1_emerg_Land  1   61.355 73.355
## mean.mud       1   63.453 75.453
## forest.dist    1   60.656 72.656
## MeanTemp       1   60.524 72.524
Bb1<-glm(Bbuf_pres~ X1_submerg + X1_emerg + X1_emerg_Land +
          mean.mud + forest.dist, 
        data, family = "binomial")
simulateResiduals(Bb1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2099373 0.2280974 0.7493404 0.9898203 0.03628088 0.06183371 0.5679061 0.06516832 0.04862612 0.166426 0.426804 0.1918813 0.5680537 0.7469708 0.9261278 0.5498223 0.5244157 0.09838691 0.9902171 0.2345293 ...
drop1(Bb1)# forest.dist 
## Single term deletions
## 
## Model:
## Bbuf_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + 
##     forest.dist
##               Df Deviance    AIC
## <none>             60.524 72.524
## X1_submerg     1   60.862 70.862
## X1_emerg       1   61.275 71.275
## X1_emerg_Land  1   61.365 71.365
## mean.mud       1   63.465 73.465
## forest.dist    1   60.662 70.662
Bb2<-glm(Bbuf_pres~ X1_submerg + X1_emerg + X1_emerg_Land +
          mean.mud, 
        data, family = "binomial")
simulateResiduals(Bb2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2082843 0.2160127 0.7396061 0.9895808 0.03668177 0.06281519 0.5498773 0.06559149 0.04780658 0.1636828 0.4395065 0.1813878 0.5530294 0.7469708 0.9234169 0.5405556 0.5403071 0.09838691 0.9892388 0.2295745 ...
drop1(Bb2)#submerg
## Single term deletions
## 
## Model:
## Bbuf_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud
##               Df Deviance    AIC
## <none>             60.662 70.662
## X1_submerg     1   60.917 68.917
## X1_emerg       1   61.305 69.305
## X1_emerg_Land  1   61.371 69.371
## mean.mud       1   63.558 71.558
Bb3<-glm(Bbuf_pres~ X1_emerg + X1_emerg_Land +
          mean.mud, 
        data, family = "binomial")
simulateResiduals(Bb3, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.1950599 0.2160127 0.7396061 0.9898203 0.03648132 0.06772263 0.5378582 0.07024637 0.04726022 0.166426 0.4344255 0.1933804 0.5718097 0.7434565 0.9234169 0.5467334 0.5371289 0.1062579 0.9894834 0.2345293 ...
drop1(Bb3)#emerg 
## Single term deletions
## 
## Model:
## Bbuf_pres ~ X1_emerg + X1_emerg_Land + mean.mud
##               Df Deviance    AIC
## <none>             60.917 68.917
## X1_emerg       1   61.307 67.307
## X1_emerg_Land  1   61.595 67.595
## mean.mud       1   63.801 69.801
summary(Bb3) # best model since without emerg and emerg_Land Quant dev is observed.
## 
## Call:
## glm(formula = Bbuf_pres ~ X1_emerg + X1_emerg_Land + mean.mud, 
##     family = "binomial", data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -0.5181     0.3099  -1.672   0.0946 .
## X1_emerg       -0.2639     0.4244  -0.622   0.5340  
## X1_emerg_Land   0.3119     0.3789   0.823   0.4105  
## mean.mud        0.6599     0.4251   1.552   0.1206  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 66.406  on 49  degrees of freedom
## Residual deviance: 60.917  on 46  degrees of freedom
## AIC: 68.917
## 
## Number of Fisher Scoring iterations: 4
tab_model(Bb3)
  Bbuf_pres
Predictors Odds Ratios CI p
(Intercept) 0.60 0.32 – 1.08 0.095
X1 emerg 0.77 0.32 – 1.76 0.534
X1 emerg Land 1.37 0.65 – 2.93 0.410
mean mud 1.93 0.91 – 4.99 0.121
Observations 50
R2 Tjur 0.105
# None of the environemental variable affect probability of presence of Bufo bufo. 

We need to extract the results from comparison between models:

# Extract AIC values
aic_values <- c(AIC(Bb), AIC(Bb1), AIC(Bb2), AIC(Bb3))
names(aic_values) <- c("Bb", "Bb1", "Bb2", "Bb3")

# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min

# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))

# Create a summary table
model_comparison <- data.frame(
  Model = names(aic_values),
  AIC = aic_values,
  Delta_AIC = delta_aic,
  Akaike_Weight = aic_weights
)
print(model_comparison)
##     Model      AIC Delta_AIC Akaike_Weight
## Bb     Bb 74.49068  5.573976    0.03746865
## Bb1   Bb1 72.52442  3.607712    0.10014676
## Bb2   Bb2 70.66152  1.744810    0.25419140
## Bb3   Bb3 68.91671  0.000000    0.60819320

Now, we need to represent it graphically:

summ(Bb3)
Observations 50
Dependent variable Bbuf_pres
Type Generalized linear model
Family binomial
Link logit
χ²(3) 5.49
Pseudo-R² (Cragg-Uhler) 0.14
Pseudo-R² (McFadden) 0.08
AIC 68.92
BIC 76.56
Est. S.E. z val. p
(Intercept) -0.52 0.31 -1.67 0.09
X1_emerg -0.26 0.42 -0.62 0.53
X1_emerg_Land 0.31 0.38 0.82 0.41
mean.mud 0.66 0.43 1.55 0.12
Standard errors: MLE
plot_summs(Bb3, colors = "#9E6969", point.size = 1, inner_ci_level = .95)
## Le chargement a nécessité le package : broom.mixed

to be able to plot all the species on the same graph, we need to extract the coeff and se of each model. using the following function.

sumbb<-summary(Bb3)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Bufo bufo"

coefBB<-dfcoef

Hyla arborea

length(data$points[data$Harbo_pres==1])
## [1] 4

PCA:

pcaH<-prcomp(data[,c(11,13:16)])
pcaH
## Standard deviations (1, .., p=5):
## [1] 0.6238767 0.5337769 0.4227561 0.3253126 0.2366494
## 
## Rotation (n x k) = (5 x 5):
##                     PC1       PC2         PC3         PC4        PC5
## Bbuf_pres    0.09524347 0.8570632  0.16128109 -0.44887571  0.1699129
## Lvulga_pres -0.16039910 0.2781172 -0.91825270  0.08550326 -0.2154626
## Pridi_pres  -0.13833213 0.3171985  0.33981549  0.40499699 -0.7750820
## Pberg_pres  -0.67017677 0.2094372  0.10692662  0.48811948  0.5072526
## Rtemp_pres   0.70493130 0.2088413 -0.06239011  0.62363178  0.2581628
str(pcaH)
## List of 5
##  $ sdev    : num [1:5] 0.624 0.534 0.423 0.325 0.237
##  $ rotation: num [1:5, 1:5] 0.0952 -0.1604 -0.1383 -0.6702 0.7049 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "Bbuf_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:5] 0.38 0.24 0.12 0.34 0.56
##   ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:50, 1:5] -0.957 -0.818 0.652 0.652 -0.148 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
PCA(data[,c(11, 13:16)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
pcaH$x[,1]#is what we use as explanatory variable. 
##  [1] -0.95650720 -0.81817507  0.65217646  0.65217646 -0.14799830 -0.14799830
##  [7] -0.14799830  0.55693299 -0.97857417 -0.30839740 -0.81817507  0.55693299
## [13]  0.65217646 -0.88333070 -0.86126373 -0.14799830 -0.81817507 -0.14799830
## [19]  0.49177736 -0.81817507 -0.05275483  0.55693299 -0.97857417 -0.86126373
## [25]  0.65217646  0.55693299  0.39653389  0.39653389 -0.05275483  0.55693299
## [31] -0.01800030 -0.81817507  0.55693299  0.51384433 -0.31673154  0.55693299
## [37]  0.55693299  0.55693299 -0.97857417  0.55693299 -0.81817507 -0.15633244
## [43]  0.55693299  0.65217646 -0.88333070  0.55693299  0.49177736  0.65217646
## [49]  0.55693299  0.49177736
Hasp<-glm(data$Harbo_pres ~ pcaH$x[,1], family = "binomial")
simulateResiduals(Hasp, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.4132625 0.3776447 0.2521923 0.9902993 0.05011171 0.1226859 0.7511986 0.09479028 0.06829512 0.2286072 0.635125 0.9374768 0.04269034 0.1214264 0.8305683 0.7722223 0.7945693 0.2459673 0.907385 0.4129037 ...
summary(Hasp)
## 
## Call:
## glm(formula = data$Harbo_pres ~ pcaH$x[, 1], family = "binomial")
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -7.807      5.461  -1.430    0.153
## pcaH$x[, 1]   10.723      9.253   1.159    0.247
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 27.877  on 49  degrees of freedom
## Residual deviance: 20.537  on 48  degrees of freedom
## AIC: 24.537
## 
## Number of Fisher Scoring iterations: 9
Ha<-glm(data$Harbo_pres~data$X1_submerg + 
          data$X1_emerg + 
          data$mean.mud + 
          data$forest.dist  + I(data$forest.dist^2)+
          data$MeanTemp + pcaH$x[,1], family = "binomial")
## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
simulateResiduals(Ha, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.4132625 0.3776447 0.391603 0.9700596 0.05011171 0.1226859 0.7511986 0.1057927 0.06829512 0.2286072 0.635125 0.3747682 0.0609862 0.1214264 0.8305683 0.7722223 0.7945693 0.2459673 0.9694284 0.4129037 ...
testZeroInflation(simulateResiduals(Ha))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1, p-value = 1
## alternative hypothesis: two.sided
summary(Ha)
## 
## Call:
## glm(formula = data$Harbo_pres ~ data$X1_submerg + data$X1_emerg + 
##     data$mean.mud + data$forest.dist + I(data$forest.dist^2) + 
##     data$MeanTemp + pcaH$x[, 1], family = "binomial")
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)
## (Intercept)             -379.84   94203.24  -0.004    0.997
## data$X1_submerg          -12.69   23524.76  -0.001    1.000
## data$X1_emerg            -55.34   27075.31  -0.002    0.998
## data$mean.mud            -34.12   22540.73  -0.002    0.999
## data$forest.dist         -32.83   15099.33  -0.002    0.998
## I(data$forest.dist^2)     29.69   12727.97   0.002    0.998
## data$MeanTemp             23.90   16517.19   0.001    0.999
## pcaH$x[, 1]              527.45  145959.86   0.004    0.997
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2.7877e+01  on 49  degrees of freedom
## Residual deviance: 4.0644e-08  on 42  degrees of freedom
## AIC: 16
## 
## Number of Fisher Scoring iterations: 25
drop1(Ha)# l'algorithme ne converge pas --> no matter the variable I try to remove. Cannot say anything about how the other sp. affect the proba of presence of H. arborea.  
## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Single term deletions
## 
## Model:
## data$Harbo_pres ~ data$X1_submerg + data$X1_emerg + data$mean.mud + 
##     data$forest.dist + I(data$forest.dist^2) + data$MeanTemp + 
##     pcaH$x[, 1]
##                       Df Deviance    AIC
## <none>                     0.0000 16.000
## data$X1_submerg        1   0.0000 14.000
## data$X1_emerg          1   0.0000 14.000
## data$mean.mud          1   0.0000 14.000
## data$forest.dist       1   9.3795 23.379
## I(data$forest.dist^2)  1   0.0000 14.000
## data$MeanTemp          1   0.0000 14.000
## pcaH$x[, 1]            1  22.5916 36.592
Ha<-glm(Harbo_pres~X1_submerg + 
          X1_emerg + 
          mean.mud + 
          forest.dist  + I(forest.dist^2)+
          MeanTemp, data, family = "binomial")
simulateResiduals(Ha, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.4000381 0.3519648 0.3665404 0.9972455 0.04971081 0.1192507 0.7271602 0.1045232 0.06610967 0.2030032 0.635125 0.8149314 0.05952253 0.1185122 0.7907011 0.7228001 0.791391 0.1839835 0.8143198 0.3253681 ...
testZeroInflation(simulateResiduals(Ha))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99896, p-value = 1
## alternative hypothesis: two.sided
summary(Ha)
## 
## Call:
## glm(formula = Harbo_pres ~ X1_submerg + X1_emerg + mean.mud + 
##     forest.dist + I(forest.dist^2) + MeanTemp, family = "binomial", 
##     data = data)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      -3.95002    1.49142  -2.648  0.00809 **
## X1_submerg        0.08516    0.71956   0.118  0.90579   
## X1_emerg         -0.84708    1.01625  -0.834  0.40454   
## mean.mud         -0.31057    0.83913  -0.370  0.71130   
## forest.dist      -1.11363    0.96179  -1.158  0.24692   
## I(forest.dist^2)  0.87463    0.92656   0.944  0.34520   
## MeanTemp          0.45575    0.99693   0.457  0.64756   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 27.877  on 49  degrees of freedom
## Residual deviance: 22.592  on 43  degrees of freedom
## AIC: 36.592
## 
## Number of Fisher Scoring iterations: 6
drop1(Ha)#submerged
## Single term deletions
## 
## Model:
## Harbo_pres ~ X1_submerg + X1_emerg + mean.mud + forest.dist + 
##     I(forest.dist^2) + MeanTemp
##                  Df Deviance    AIC
## <none>                22.592 36.592
## X1_submerg        1   22.606 34.606
## X1_emerg          1   23.252 35.252
## mean.mud          1   22.747 34.747
## forest.dist       1   24.165 36.165
## I(forest.dist^2)  1   23.448 35.448
## MeanTemp          1   22.799 34.799
Ha1<-glm(Harbo_pres~ X1_emerg + 
          mean.mud + 
          forest.dist  + I(forest.dist^2)+
          MeanTemp, data, family = "binomial")
simulateResiduals(Ha1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.4016911 0.3519648 0.3696732 0.9973652 0.04971081 0.1192507 0.730165 0.1041 0.06610967 0.2030032 0.635125 0.8099295 0.05952253 0.1185122 0.7907011 0.7228001 0.791391 0.1829996 0.8181976 0.3253681 ...
drop1(Ha1)#mean.mud
## Single term deletions
## 
## Model:
## Harbo_pres ~ X1_emerg + mean.mud + forest.dist + I(forest.dist^2) + 
##     MeanTemp
##                  Df Deviance    AIC
## <none>                22.606 34.606
## X1_emerg          1   23.354 33.354
## mean.mud          1   22.765 32.765
## forest.dist       1   24.265 34.265
## I(forest.dist^2)  1   23.470 33.470
## MeanTemp          1   22.921 32.921
Ha2<-glm(Harbo_pres~ X1_emerg + 
          forest.dist  + I(forest.dist^2)+
          MeanTemp, data, family = "binomial")
simulateResiduals(Ha2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.4000381 0.3489437 0.3681068 0.9976048 0.04971081 0.1192507 0.730165 0.1049464 0.06583649 0.2048321 0.635125 0.7699147 0.05927859 0.1189979 0.7740897 0.7320667 0.791391 0.1829996 0.8414638 0.3187617 ...
drop1(Ha2)#temp
## Single term deletions
## 
## Model:
## Harbo_pres ~ X1_emerg + forest.dist + I(forest.dist^2) + MeanTemp
##                  Df Deviance    AIC
## <none>                22.765 32.765
## X1_emerg          1   23.368 31.368
## forest.dist       1   24.448 32.448
## I(forest.dist^2)  1   23.564 31.564
## MeanTemp          1   23.047 31.047
Ha3<-glm(Harbo_pres~ X1_emerg + 
          forest.dist  + I(forest.dist^2), 
         data, family = "binomial")
simulateResiduals(Ha3, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.4000381 0.3459225 0.364974 0.9976048 0.04971081 0.1192507 0.7271602 0.1049464 0.06556331 0.2066609 0.635125 0.8674508 0.05854676 0.1199693 0.7740897 0.7228001 0.791391 0.1623384 0.8065644 0.3270198 ...
drop1(Ha3)# forest
## Single term deletions
## 
## Model:
## Harbo_pres ~ X1_emerg + forest.dist + I(forest.dist^2)
##                  Df Deviance    AIC
## <none>                23.047 31.047
## X1_emerg          1   27.061 33.061
## forest.dist       1   24.521 30.521
## I(forest.dist^2)  1   23.813 29.813
summary(Ha3)# best model since without forest dist Quant dev is observed.
## 
## Call:
## glm(formula = Harbo_pres ~ X1_emerg + forest.dist + I(forest.dist^2), 
##     family = "binomial", data = data)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       -3.8411     1.4256  -2.694  0.00705 **
## X1_emerg          -1.1193     0.6210  -1.802  0.07148 . 
## forest.dist       -0.9895     0.8586  -1.153  0.24909   
## I(forest.dist^2)   0.8155     0.8940   0.912  0.36167   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 27.877  on 49  degrees of freedom
## Residual deviance: 23.047  on 46  degrees of freedom
## AIC: 31.047
## 
## Number of Fisher Scoring iterations: 6
#hyla arbo tends to be affected by proportion of emerged vegetation. 

We need to extract the results from comparison between models:

# Extract AIC values
aic_values <- c(AIC(Ha), AIC(Ha1), AIC(Ha2), AIC(Ha3))
names(aic_values) <- c("Ha", "Ha1", "Ha2", "Ha3")

# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min

# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))

# Create a summary table
model_comparison <- data.frame(
  Model = names(aic_values),
  AIC = aic_values,
  Delta_AIC = delta_aic,
  Akaike_Weight = aic_weights
)
print(model_comparison)
##     Model      AIC Delta_AIC Akaike_Weight
## Ha     Ha 36.59158  5.544871    0.03777518
## Ha1   Ha1 34.60555  3.558839    0.10196894
## Ha2   Ha2 32.76496  1.718247    0.25594550
## Ha3   Ha3 31.04671  0.000000    0.60431037

Graphic representaiton

summ(Ha3)
Observations 50
Dependent variable Harbo_pres
Type Generalized linear model
Family binomial
Link logit
χ²(3) 4.83
Pseudo-R² (Cragg-Uhler) 0.22
Pseudo-R² (McFadden) 0.17
AIC 31.05
BIC 38.69
Est. S.E. z val. p
(Intercept) -3.84 1.43 -2.69 0.01
X1_emerg -1.12 0.62 -1.80 0.07
forest.dist -0.99 0.86 -1.15 0.25
I(forest.dist^2) 0.82 0.89 0.91 0.36
Standard errors: MLE
plot_summs(Ha3, inner_ci_level = .95, colors = "black", plot.distributions = T)

sumbb<-summary(Ha3)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Hyla arborea"

coefHA<-dfcoef

coef_final<- full_join(coefBB, coefHA)
## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`

Lissotriton vulgaris

length(data$points[data$Lvulga_pres==1])
## [1] 12

PCA

pcaL<-prcomp(data[,c(11:12,14:16)])
pcaL
## Standard deviations (1, .., p=5):
## [1] 0.6255488 0.5287522 0.3263076 0.2665902 0.2350240
## 
## Rotation (n x k) = (5 x 5):
##                   PC1       PC2           PC3        PC4         PC5
## Bbuf_pres   0.1688900 0.8570597  0.4461803382  0.1936511 -0.01863562
## Harbo_pres  0.1505903 0.1244790  0.0001832353 -0.7441496 -0.63880266
## Pridi_pres -0.1177651 0.3855602 -0.4560878722 -0.4930396  0.62158714
## Pberg_pres -0.6470553 0.3001270 -0.4741697480  0.2872603 -0.42882144
## Rtemp_pres  0.7185051 0.1059285 -0.6066877081  0.2883298 -0.14603246
str(pcaL)
## List of 5
##  $ sdev    : num [1:5] 0.626 0.529 0.326 0.267 0.235
##  $ rotation: num [1:5, 1:5] 0.169 0.151 -0.118 -0.647 0.719 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Pridi_pres" "Pberg_pres" ...
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:5] 0.38 0.08 0.12 0.34 0.56
##   ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Pridi_pres" "Pberg_pres" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:50, 1:5] -1.009 -0.892 0.643 0.794 -0.244 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
PCA(data[,c(11:12, 14:16)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
pcaL$x[,1]
##  [1] -1.009278032 -0.891512941  0.642937410  0.793527677 -0.244457678
##  [6] -0.244457678 -0.244457678  0.474047454 -0.891512941 -0.244457678
## [11] -0.891512941  0.624637721  0.642937410 -0.722622985 -0.840388077
## [16] -0.244457678 -0.891512941 -0.244457678  0.642937410 -0.891512941
## [21] -0.075567723  0.474047454 -0.891512941 -0.840388077  0.642937410
## [26]  0.474047454  0.474047454  0.474047454 -0.075567723  0.474047454
## [31] -0.004117853 -0.891512941  0.474047454  0.675762585 -0.121882945
## [36]  0.474047454  0.474047454  0.474047454 -0.891512941  0.474047454
## [41] -0.891512941 -0.121882945  0.474047454  0.793527677 -0.722622985
## [46]  0.474047454  0.642937410  0.642937410  0.474047454  0.642937410
Lvsp<-glm(data$Lvulga_pres ~ pcaL$x[,1], family = "binomial")
simulateResiduals(Lvsp, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2578758 0.2583089 0.2991847 0.7605268 0.03928758 0.09471353 0.5528821 0.0871732 0.7167617 0.8210369 0.4293445 0.3133062 0.05025263 0.761028 0.5946869 0.5621778 0.5561985 0.1859513 0.9953531 0.2906842 ...
summary(Lvsp)
## 
## Call:
## glm(formula = data$Lvulga_pres ~ pcaL$x[, 1], family = "binomial")
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1730     0.3369  -3.481 0.000499 ***
## pcaL$x[, 1]  -0.4533     0.5302  -0.855 0.392574    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.108  on 49  degrees of freedom
## Residual deviance: 54.376  on 48  degrees of freedom
## AIC: 58.376
## 
## Number of Fisher Scoring iterations: 4
Lv<-glm(data$Lvulga_pres~data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +
          data$mean.mud + data$forest.dist + data$MeanTemp + pcaL$x[,1], family = "binomial")
simulateResiduals(Lv, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.332263 0.2462243 0.3696732 0.8885746 0.04730545 0.1035469 0.5889397 0.09394394 0.6012303 0.7963523 0.579234 0.1394138 0.05830281 0.8137424 0.3056492 0.6486667 0.7596083 0.1938222 0.9953531 0.115613 ...
testZeroInflation(simulateResiduals(Lv))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99435, p-value = 1
## alternative hypothesis: two.sided
summary(Lv)
## 
## Call:
## glm(formula = data$Lvulga_pres ~ data$X1_submerg + data$X1_emerg + 
##     data$X1_emerg_Land + data$mean.mud + data$forest.dist + data$MeanTemp + 
##     pcaL$x[, 1], family = "binomial")
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -1.44773    0.41601  -3.480 0.000501 ***
## data$X1_submerg     0.20681    0.49568   0.417 0.676509    
## data$X1_emerg       1.28867    0.94783   1.360 0.173958    
## data$X1_emerg_Land -0.80793    0.51049  -1.583 0.113501    
## data$mean.mud       0.08648    0.47946   0.180 0.856855    
## data$forest.dist   -0.86469    0.49556  -1.745 0.081009 .  
## data$MeanTemp       1.32612    0.78721   1.685 0.092071 .  
## pcaL$x[, 1]        -0.74118    0.65698  -1.128 0.259257    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.108  on 49  degrees of freedom
## Residual deviance: 46.540  on 42  degrees of freedom
## AIC: 62.54
## 
## Number of Fisher Scoring iterations: 5
drop1(Lv)# mean.mud
## Single term deletions
## 
## Model:
## data$Lvulga_pres ~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land + 
##     data$mean.mud + data$forest.dist + data$MeanTemp + pcaL$x[, 
##     1]
##                    Df Deviance    AIC
## <none>                  46.540 62.540
## data$X1_submerg     1   46.715 60.715
## data$X1_emerg       1   48.731 62.731
## data$X1_emerg_Land  1   49.464 63.464
## data$mean.mud       1   46.572 60.572
## data$forest.dist    1   50.117 64.117
## data$MeanTemp       1   49.813 63.813
## pcaL$x[, 1]         1   47.870 61.870
Lv1<-glm(data$Lvulga_pres~data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +data$forest.dist + data$MeanTemp + pcaL$x[,1], family = "binomial")
simulateResiduals(Lv1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3471405 0.2492455 0.3696732 0.8885746 0.04730545 0.1040377 0.5889397 0.09267443 0.6086839 0.7932667 0.579234 0.1469091 0.05854676 0.8102281 0.315616 0.6455778 0.7596083 0.1918545 0.9952308 0.1172647 ...
drop1(Lv1)#submerg
## Single term deletions
## 
## Model:
## data$Lvulga_pres ~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land + 
##     data$forest.dist + data$MeanTemp + pcaL$x[, 1]
##                    Df Deviance    AIC
## <none>                  46.572 60.572
## data$X1_submerg     1   46.742 58.742
## data$X1_emerg       1   48.866 60.866
## data$X1_emerg_Land  1   49.548 61.548
## data$forest.dist    1   50.119 62.119
## data$MeanTemp       1   49.825 61.825
## pcaL$x[, 1]         1   47.932 59.932
Lv2<-glm(data$Lvulga_pres~data$X1_emerg + data$X1_emerg_Land +data$forest.dist + data$MeanTemp + pcaL$x[,1], family = "binomial")
simulateResiduals(Lv2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3603649 0.2492455 0.3712397 0.8846944 0.04730545 0.1006025 0.6009588 0.08971223 0.627318 0.7901812 0.579234 0.1379147 0.05854676 0.8102281 0.315616 0.6424889 0.7500734 0.1898867 0.9952308 0.107355 ...
drop1(Lv2)# pca
## Single term deletions
## 
## Model:
## data$Lvulga_pres ~ data$X1_emerg + data$X1_emerg_Land + data$forest.dist + 
##     data$MeanTemp + pcaL$x[, 1]
##                    Df Deviance    AIC
## <none>                  46.742 58.742
## data$X1_emerg       1   48.885 58.885
## data$X1_emerg_Land  1   49.623 59.623
## data$forest.dist    1   50.122 60.122
## data$MeanTemp       1   50.406 60.406
## pcaL$x[, 1]         1   48.143 58.143
Lv3<-glm(data$Lvulga_pres~data$X1_emerg + data$X1_emerg_Land +data$forest.dist + data$MeanTemp, family = "binomial")
simulateResiduals(Lv3, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3752423 0.2854994 0.3524427 0.8148501 0.047105 0.1055099 0.5829301 0.08548052 0.7055813 0.7963523 0.604639 0.09144343 0.05537547 0.845371 0.4318955 0.6424889 0.7754997 0.1987416 0.9919291 0.1701163 ...
drop1(Lv3)# best model --> sp. does not affect the probability of presence of L. vulgaris. 
## Single term deletions
## 
## Model:
## data$Lvulga_pres ~ data$X1_emerg + data$X1_emerg_Land + data$forest.dist + 
##     data$MeanTemp
##                    Df Deviance    AIC
## <none>                  48.143 58.143
## data$X1_emerg       1   50.632 58.632
## data$X1_emerg_Land  1   51.149 59.149
## data$forest.dist    1   52.117 60.117
## data$MeanTemp       1   51.303 59.303
summary(Lv3)
## 
## Call:
## glm(formula = data$Lvulga_pres ~ data$X1_emerg + data$X1_emerg_Land + 
##     data$forest.dist + data$MeanTemp, family = "binomial")
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.3826     0.3974  -3.479 0.000503 ***
## data$X1_emerg        1.1799     0.8163   1.445 0.148337    
## data$X1_emerg_Land  -0.8353     0.5218  -1.601 0.109414    
## data$forest.dist    -0.8669     0.4848  -1.788 0.073716 .  
## data$MeanTemp        1.2223     0.7487   1.632 0.102582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.108  on 49  degrees of freedom
## Residual deviance: 48.143  on 45  degrees of freedom
## AIC: 58.143
## 
## Number of Fisher Scoring iterations: 5
Lv<-glm(Lvulga_pres~X1_submerg + X1_emerg + X1_emerg_Land +
          mean.mud + forest.dist + MeanTemp, 
        data, family = "binomial")
simulateResiduals(Lv, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3603649 0.2854994 0.3493099 0.8187303 0.04690456 0.1074729 0.5528821 0.0901354 0.649679 0.8086946 0.604639 0.0929425 0.05464364 0.8418568 0.4119619 0.6486667 0.7754997 0.2075964 0.9930297 0.1783744 ...
testZeroInflation(simulateResiduals(Lv))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99643, p-value = 1
## alternative hypothesis: two.sided
drop1(Lv)#mean.mud
## Single term deletions
## 
## Model:
## Lvulga_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + 
##     forest.dist + MeanTemp
##               Df Deviance    AIC
## <none>             47.870 61.870
## X1_submerg     1   48.082 60.082
## X1_emerg       1   50.500 62.500
## X1_emerg_Land  1   50.905 62.905
## mean.mud       1   47.932 59.932
## forest.dist    1   52.092 64.092
## MeanTemp       1   50.630 62.630
Lv1<-glm(Lvulga_pres~X1_submerg + X1_emerg + X1_emerg_Land +
          forest.dist + MeanTemp, 
        data, family = "binomial")
simulateResiduals(Lv, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3603649 0.2854994 0.3493099 0.8187303 0.04690456 0.1074729 0.5528821 0.0901354 0.649679 0.8086946 0.604639 0.0929425 0.05464364 0.8418568 0.4119619 0.6486667 0.7754997 0.2075964 0.9930297 0.1783744 ...
drop1(Lv1)#submerg
## Single term deletions
## 
## Model:
## Lvulga_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + forest.dist + 
##     MeanTemp
##               Df Deviance    AIC
## <none>             47.932 59.932
## X1_submerg     1   48.143 58.143
## X1_emerg       1   50.622 60.622
## X1_emerg_Land  1   51.014 61.014
## forest.dist    1   52.095 62.095
## MeanTemp       1   50.656 60.656
Lv2<-glm(Lvulga_pres~X1_emerg + X1_emerg_Land +
          forest.dist + MeanTemp, 
        data, family = "binomial")
simulateResiduals(Lv2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3752423 0.2854994 0.3524427 0.8148501 0.047105 0.1055099 0.5829301 0.08548052 0.7055813 0.7963523 0.604639 0.09144343 0.05537547 0.845371 0.4318955 0.6424889 0.7754997 0.1987416 0.9919291 0.1701163 ...
drop1(Lv2)# best model 
## Single term deletions
## 
## Model:
## Lvulga_pres ~ X1_emerg + X1_emerg_Land + forest.dist + MeanTemp
##               Df Deviance    AIC
## <none>             48.143 58.143
## X1_emerg       1   50.632 58.632
## X1_emerg_Land  1   51.149 59.149
## forest.dist    1   52.117 60.117
## MeanTemp       1   51.303 59.303
summary(Lv2)
## 
## Call:
## glm(formula = Lvulga_pres ~ X1_emerg + X1_emerg_Land + forest.dist + 
##     MeanTemp, family = "binomial", data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.3826     0.3974  -3.479 0.000503 ***
## X1_emerg        1.1799     0.8163   1.445 0.148337    
## X1_emerg_Land  -0.8353     0.5218  -1.601 0.109414    
## forest.dist    -0.8669     0.4848  -1.788 0.073716 .  
## MeanTemp        1.2223     0.7487   1.632 0.102582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.108  on 49  degrees of freedom
## Residual deviance: 48.143  on 45  degrees of freedom
## AIC: 58.143
## 
## Number of Fisher Scoring iterations: 5
#L. vulgaris tends to be affect by the dist to the wintering habitats. 

We need to extract the results from comparison between models:

# Extract AIC values
aic_values <- c(AIC(Lv), AIC(Lv1), AIC(Lv2))
names(aic_values) <- c("Lv", "Lv1", "Lv2")

# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min

# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))

# Create a summary table
model_comparison <- data.frame(
  Model = names(aic_values),
  AIC = aic_values,
  Delta_AIC = delta_aic,
  Akaike_Weight = aic_weights
)
print(model_comparison)
##     Model      AIC Delta_AIC Akaike_Weight
## Lv     Lv 61.86963  3.726885    0.09919995
## Lv1   Lv1 59.93204  1.789299    0.26136830
## Lv2   Lv2 58.14274  0.000000    0.63943175

Graphical representation:

summ(Lv2)
Observations 50
Dependent variable Lvulga_pres
Type Generalized linear model
Family binomial
Link logit
χ²(4) 6.97
Pseudo-R² (Cragg-Uhler) 0.19
Pseudo-R² (McFadden) 0.13
AIC 58.14
BIC 67.70
Est. S.E. z val. p
(Intercept) -1.38 0.40 -3.48 0.00
X1_emerg 1.18 0.82 1.45 0.15
X1_emerg_Land -0.84 0.52 -1.60 0.11
forest.dist -0.87 0.48 -1.79 0.07
MeanTemp 1.22 0.75 1.63 0.10
Standard errors: MLE
plot_summs(Lv2, inner_ci_level = .95, colors = "black", plot.distributions = T)

sumbb<-summary(Lv2)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Lissotriton vulgaris"

coefLv<-dfcoef

coef_final<- full_join(coef_final, coefLv)
## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`

Pelophylax ridibundus

length(data$points[data$Pridi_pres==1])
## [1] 6

PCA

pcaPr<-prcomp(data[,c(11:13,15:16)])
pcaPr
## Standard deviations (1, .., p=5):
## [1] 0.6269371 0.5159911 0.4089307 0.3066318 0.2479851
## 
## Rotation (n x k) = (5 x 5):
##                    PC1        PC2        PC3         PC4         PC5
## Bbuf_pres    0.1814447 0.85962133  0.3745276 -0.23949634  0.17464118
## Harbo_pres   0.1633600 0.07269459  0.1857509 -0.02906505 -0.96575403
## Lvulga_pres -0.1606287 0.41296865 -0.8734737 -0.12234836 -0.16040537
## Pberg_pres  -0.6342752 0.27393148  0.1621125  0.70036667 -0.07656764
## Rtemp_pres   0.7157407 0.10092124 -0.1897073  0.66054078  0.07229882
str(pcaPr)
## List of 5
##  $ sdev    : num [1:5] 0.627 0.516 0.409 0.307 0.248
##  $ rotation: num [1:5, 1:5] 0.181 0.163 -0.161 -0.634 0.716 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pberg_pres" ...
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:5] 0.38 0.08 0.24 0.34 0.56
##   ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pberg_pres" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:50, 1:5] -0.863 -0.863 0.669 0.832 -0.229 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
PCA(data[,c(11:13, 15:16)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
pcaPr$x[,1]
##  [1] -0.86290331 -0.86290331  0.66855725  0.83191725 -0.22862814 -0.22862814
##  [7] -0.22862814  0.48711258 -1.02353198 -0.38925680 -0.86290331  0.65047258
## [13]  0.66855725 -0.84208731 -0.68145865 -0.22862814 -0.86290331 -0.22862814
## [19]  0.50792859 -0.86290331 -0.04718347  0.48711258 -1.02353198 -0.68145865
## [25]  0.66855725  0.48711258  0.32648392  0.32648392 -0.04718347  0.48711258
## [31]  0.03428208 -0.86290331  0.48711258  0.83191725 -0.12634658  0.48711258
## [37]  0.48711258  0.48711258 -1.02353198  0.48711258 -0.86290331  0.03428208
## [43]  0.48711258  0.83191725 -0.84208731  0.48711258  0.50792859  0.66855725
## [49]  0.48711258  0.50792859
Prsp<-glm(data$Pridi_pres~ pcaPr$x[,1], family = "binomial")
simulateResiduals(Prsp, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.859183 0.3021157 0.3602748 0.8846944 0.04550143 0.1133618 0.6400212 0.1011378 0.05545564 0.1956878 0.5081 0.3552802 0.05781492 0.09714111 0.9722132 0.6888223 0.6547251 0.2095641 0.8957518 0.3270198 ...
summary(Prsp)
## 
## Call:
## glm(formula = data$Pridi_pres ~ pcaPr$x[, 1], family = "binomial")
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.0678     0.4641  -4.455 8.38e-06 ***
## pcaPr$x[, 1]  -0.7130     0.6982  -1.021    0.307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 36.692  on 49  degrees of freedom
## Residual deviance: 35.632  on 48  degrees of freedom
## AIC: 39.632
## 
## Number of Fisher Scoring iterations: 5
Pr<-glm(Pridi_pres~X1_emerg + 
          X1_emerg_Land + 
          mean.mud + I(mean.mud^2) + 
          forest.dist  +   
          MeanTemp + pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.624488 0.2583089 0.3712397 0.9467782 0.04770634 0.1109081 0.7211506 0.1053696 0.04807976 0.2194629 0.5817745 0.3417886 0.05561942 0.1131694 0.9390046 0.7567778 0.7500734 0.1889029 0.9422844 0.2807745 ...
testZeroInflation(simulateResiduals(Pr))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99233, p-value = 1
## alternative hypothesis: two.sided
drop1(Pr)#emerged vegetation
## Single term deletions
## 
## Model:
## Pridi_pres ~ X1_emerg + X1_emerg_Land + mean.mud + I(mean.mud^2) + 
##     forest.dist + MeanTemp + pcaPr$x[, 1]
##               Df Deviance    AIC
## <none>             28.549 44.549
## X1_emerg       1   28.550 42.550
## X1_emerg_Land  1   28.812 42.812
## mean.mud       1   33.434 47.434
## I(mean.mud^2)  1   30.405 44.405
## forest.dist    1   28.557 42.557
## MeanTemp       1   28.569 42.569
## pcaPr$x[, 1]   1   29.897 43.897
Pr1<-glm(Pridi_pres~X1_emerg_Land + 
          mean.mud + I(mean.mud^2) + 
          forest.dist  +   
          MeanTemp + pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.624488 0.2583089 0.3712397 0.9467782 0.04770634 0.1109081 0.7211506 0.1053696 0.04807976 0.2194629 0.584315 0.3417886 0.05561942 0.1131694 0.9383269 0.7567778 0.7500734 0.1889029 0.9422844 0.2807745 ...
drop1(Pr1)#forest dist
## Single term deletions
## 
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + forest.dist + 
##     MeanTemp + pcaPr$x[, 1]
##               Df Deviance    AIC
## <none>             28.550 42.550
## X1_emerg_Land  1   28.845 40.845
## mean.mud       1   33.649 45.649
## I(mean.mud^2)  1   30.406 42.406
## forest.dist    1   28.558 40.558
## MeanTemp       1   28.581 40.581
## pcaPr$x[, 1]   1   29.928 41.928
Pr2<-glm(Pridi_pres~X1_emerg_Land + 
          mean.mud + I(mean.mud^2) +    
          MeanTemp + pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.6291819 0.2583089 0.3696732 0.9467782 0.04770634 0.1099266 0.7211506 0.1053696 0.04807976 0.2194629 0.584315 0.3402895 0.05537547 0.1131694 0.9376492 0.7567778 0.7532517 0.1889029 0.9422844 0.2758197 ...
drop1(Pr2)#Mean temp
## Single term deletions
## 
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + MeanTemp + 
##     pcaPr$x[, 1]
##               Df Deviance    AIC
## <none>             28.558 40.558
## X1_emerg_Land  1   28.847 38.847
## mean.mud       1   34.353 44.353
## I(mean.mud^2)  1   30.742 40.742
## MeanTemp       1   28.582 38.582
## pcaPr$x[, 1]   1   29.933 39.933
Pr3<-glm(Pridi_pres~X1_emerg_Land + 
          mean.mud + I(mean.mud^2) + 
           pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr3, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.6338758 0.2583089 0.3712397 0.9467782 0.04790679 0.1069821 0.7211506 0.1053696 0.04862612 0.2185485 0.579234 0.3402895 0.05610731 0.112198 0.9383269 0.7567778 0.7405386 0.1889029 0.9422844 0.2725165 ...
drop1(Pr3)#emerg land
## Single term deletions
## 
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + pcaPr$x[, 
##     1]
##               Df Deviance    AIC
## <none>             28.582 38.582
## X1_emerg_Land  1   29.014 37.014
## mean.mud       1   34.569 42.569
## I(mean.mud^2)  1   30.743 38.743
## pcaPr$x[, 1]   1   30.320 38.320
summary(Pr3)
## 
## Call:
## glm(formula = Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + 
##     pcaPr$x[, 1], family = "binomial", data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -2.2974     0.6249  -3.677 0.000236 ***
## X1_emerg_Land   0.3788     0.5774   0.656 0.511875    
## mean.mud        2.0199     0.9811   2.059 0.039510 *  
## I(mean.mud^2)  -0.3983     0.3000  -1.327 0.184346    
## pcaPr$x[, 1]   -1.0554     0.8243  -1.280 0.200461    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 36.692  on 49  degrees of freedom
## Residual deviance: 28.582  on 45  degrees of freedom
## AIC: 38.582
## 
## Number of Fisher Scoring iterations: 6
Pr4<-glm(Pridi_pres~X1_emerg_Land + mean.mud + I(mean.mud^2), data, family = "binomial")
simulateResiduals(Pr4, plot=T)#model not converging. 

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.7066312 0.308158 0.3289465 0.9196165 0.04770634 0.1045284 0.7211506 0.1045232 0.05955334 0.2194629 0.604639 0.2923192 0.04756924 0.1146265 0.9647582 0.7567778 0.7659648 0.1918545 0.9151404 0.3270198 ...
drop1(Pr4)
## Single term deletions
## 
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2)
##               Df Deviance    AIC
## <none>             30.320 38.320
## X1_emerg_Land  1   31.122 37.122
## mean.mud       1   36.139 42.139
## I(mean.mud^2)  1   32.955 38.955

We need to extract the results from comparison between models:

# Extract AIC values
aic_values <- c(AIC(Pr), AIC(Pr1), AIC(Pr2), AIC(Pr3))
names(aic_values) <- c("Pr", "Pr1", "Pr2", "Pr3")

# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min

# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))

# Create a summary table
model_comparison <- data.frame(
  Model = names(aic_values),
  AIC = aic_values,
  Delta_AIC = delta_aic,
  Akaike_Weight = aic_weights
)
print(model_comparison)
##     Model      AIC Delta_AIC Akaike_Weight
## Pr     Pr 44.54933  5.967068    0.03243256
## Pr1   Pr1 42.55035  3.968080    0.08811621
## Pr2   Pr2 40.55754  1.975276    0.23866451
## Pr3   Pr3 38.58227  0.000000    0.64078671
summ(Pr3)
Observations 50
Dependent variable Pridi_pres
Type Generalized linear model
Family binomial
Link logit
χ²(4) 8.11
Pseudo-R² (Cragg-Uhler) 0.29
Pseudo-R² (McFadden) 0.22
AIC 38.58
BIC 48.14
Est. S.E. z val. p
(Intercept) -2.30 0.62 -3.68 0.00
X1_emerg_Land 0.38 0.58 0.66 0.51
mean.mud 2.02 0.98 2.06 0.04
I(mean.mud^2) -0.40 0.30 -1.33 0.18
pcaPr$x[, 1] -1.06 0.82 -1.28 0.20
Standard errors: MLE
plot_summs(Pr3, inner_ci_level = .95, colors = "Black", plot.distributions = T)

sumbb<-summary(Pr3)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Pelophylax ridibundus"

coefPr<-dfcoef

coef_final<- full_join(coef_final, coefPr)
## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`

Pelophylax bergeri

length(data$points[data$Pberg_pres==1])
## [1] 17

PCA

pcaPb<-prcomp(data[,c(11:14,16)])
pcaPb
## Standard deviations (1, .., p=5):
## [1] 0.5568887 0.4812606 0.4214406 0.2822512 0.2485994
## 
## Rotation (n x k) = (5 x 5):
##                    PC1         PC2        PC3         PC4         PC5
## Bbuf_pres   0.70078214  0.52160158  0.2592054  0.38696777  0.14108400
## Harbo_pres  0.19546022 -0.09325664  0.1415840  0.02890934 -0.96551371
## Lvulga_pres 0.04402947  0.53328906 -0.8038585 -0.19942801 -0.16644566
## Pridi_pres  0.17650508  0.22390988  0.3688108 -0.88371245  0.04172784
## Rtemp_pres  0.66151876 -0.62024361 -0.3613262 -0.16941362  0.13576896
str(pcaPb)
## List of 5
##  $ sdev    : num [1:5] 0.557 0.481 0.421 0.282 0.249
##  $ rotation: num [1:5, 1:5] 0.701 0.195 0.044 0.177 0.662 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:5] 0.38 0.08 0.24 0.12 0.56
##   ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:50, 1:5] -0.508 -0.684 0.678 0.874 -0.684 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
PCA(data[,c(11:14,16)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
pcaPb$x[,1]
##  [1] -0.50762714 -0.68413222  0.67816868  0.87362889 -0.68413222 -0.68413222
##  [7] -0.68413222 -0.02261346 -0.64010274 -0.64010274 -0.68413222  0.17284676
## [13]  0.67816868  0.06067939  0.19315500 -0.68413222 -0.68413222 -0.68413222
## [19]  0.72219815 -0.68413222  0.01664992 -0.02261346 -0.64010274  0.19315500
## [25]  0.67816868 -0.02261346  0.02141601  0.02141601  0.01664992 -0.02261346
## [31]  0.67816868 -0.68413222 -0.02261346  1.05013397  0.89870323 -0.02261346
## [37] -0.02261346 -0.02261346 -0.64010274 -0.02261346 -0.68413222  0.85467376
## [43] -0.02261346  0.87362889  0.06067939 -0.02261346  0.72219815  0.67816868
## [49] -0.02261346  0.72219815
Pbsp<-glm(Pberg_pres ~ pcaPb$x[,1], data, family = "binomial")
simulateResiduals(Pbsp, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.7042843 0.6838435 0.3023175 0.7954489 0.02245004 0.05741701 0.3846137 0.07363174 0.5304207 0.115218 0.8117245 0.2713321 0.05049658 0.7504851 0.9566255 0.3459556 0.8857805 0.1151127 0.8492193 0.6735745 ...
summary(Pbsp)
## 
## Call:
## glm(formula = Pberg_pres ~ pcaPb$x[, 1], family = "binomial", 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -0.7280     0.3181  -2.289   0.0221 *
## pcaPb$x[, 1]  -1.1420     0.6046  -1.889   0.0589 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 64.104  on 49  degrees of freedom
## Residual deviance: 60.150  on 48  degrees of freedom
## AIC: 64.15
## 
## Number of Fisher Scoring iterations: 4
Pb<-glm(Pberg_pres~X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + forest.dist +  MeanTemp + pcaPb$x[,1], 
        data, family = "binomial")
simulateResiduals(Pb, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.5728551 0.6863329 0.2913526 0.8769339 0.01743887 0.07459304 0.291465 0.0930976 0.3366261 0.147223 0.7854535 0.2953173 0.05074052 0.6837135 0.9471373 0.4973111 0.8701678 0.2115319 0.9384067 0.7956905 ...
testZeroInflation(simulateResiduals(Pb))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99338, p-value = 1
## alternative hypothesis: two.sided
drop1(Pb) # MeanTemp
## Single term deletions
## 
## Model:
## Pberg_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + 
##     forest.dist + MeanTemp + pcaPb$x[, 1]
##               Df Deviance    AIC
## <none>             53.358 69.358
## X1_submerg     1   53.579 67.579
## X1_emerg       1   56.597 70.597
## X1_emerg_Land  1   53.556 67.556
## mean.mud       1   55.692 69.692
## forest.dist    1   53.571 67.571
## MeanTemp       1   53.440 67.440
## pcaPb$x[, 1]   1   56.650 70.650
Pb1<-glm(Pberg_pres~X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + forest.dist  + pcaPb$x[,1], 
        data, family = "binomial")
simulateResiduals(Pb1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.5681612 0.6938012 0.2897862 0.8769339 0.01663709 0.07802825 0.2884602 0.09563662 0.3291725 0.1499663 0.7796155 0.308809 0.04976474 0.7083136 0.9498482 0.4880445 0.8709895 0.2085802 0.934529 0.7956905 ...
drop1(Pb1) # emerg Land
## Single term deletions
## 
## Model:
## Pberg_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + 
##     forest.dist + pcaPb$x[, 1]
##               Df Deviance    AIC
## <none>             53.440 67.440
## X1_submerg     1   53.665 65.665
## X1_emerg       1   58.010 70.010
## X1_emerg_Land  1   53.602 65.602
## mean.mud       1   55.791 67.791
## forest.dist    1   53.837 65.837
## pcaPb$x[, 1]   1   56.694 68.694
Pb2<-glm(Pberg_pres~X1_submerg + X1_emerg + mean.mud + forest.dist  + pcaPb$x[,1], 
        data, family = "binomial")
simulateResiduals(Pb2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.5399978 0.6888223 0.2788213 0.8536525 0.01543441 0.06870412 0.2944698 0.09225126 0.3440798 0.1490519 0.7591825 0.3133062 0.04854502 0.6766849 0.9525591 0.4787778 0.8594854 0.2085802 0.9384067 0.8074324 ...
drop1(Pb2) # submerged
## Single term deletions
## 
## Model:
## Pberg_pres ~ X1_submerg + X1_emerg + mean.mud + forest.dist + 
##     pcaPb$x[, 1]
##              Df Deviance    AIC
## <none>            53.602 65.602
## X1_submerg    1   53.816 63.816
## X1_emerg      1   58.226 68.226
## mean.mud      1   56.040 66.040
## forest.dist   1   54.178 64.178
## pcaPb$x[, 1]  1   56.846 66.846
Pb3<-glm(Pberg_pres~X1_emerg + mean.mud + forest.dist  + pcaPb$x[,1], 
        data, family = "binomial")
simulateResiduals(Pb3, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.575202 0.6863329 0.2819542 0.8575327 0.01463262 0.0662504 0.3094938 0.08928906 0.3440798 0.1453942 0.757723 0.308809 0.04854502 0.6837135 0.9545923 0.4756889 0.8594854 0.2046448 0.9384067 0.7956905 ...
drop1(Pb3) # forst.dist
## Single term deletions
## 
## Model:
## Pberg_pres ~ X1_emerg + mean.mud + forest.dist + pcaPb$x[, 1]
##              Df Deviance    AIC
## <none>            53.816 63.816
## X1_emerg      1   59.208 67.208
## mean.mud      1   56.360 64.360
## forest.dist   1   54.721 62.721
## pcaPb$x[, 1]  1   57.417 65.417
Pb4<-glm(Pberg_pres~X1_emerg + mean.mud + pcaPb$x[,1], 
        data, family = "binomial")
simulateResiduals(Pb4, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.5540795 0.6788646 0.2474931 0.8032094 0.01603575 0.05987073 0.2644219 0.08759638 0.3366261 0.1316778 0.7708585 0.2788275 0.04317823 0.6485706 0.9423932 0.4417111 0.863594 0.1967738 0.907385 0.7839486 ...
drop1(Pb4)# bestmodel
## Single term deletions
## 
## Model:
## Pberg_pres ~ X1_emerg + mean.mud + pcaPb$x[, 1]
##              Df Deviance    AIC
## <none>            54.721 62.721
## X1_emerg      1   59.997 65.997
## mean.mud      1   57.065 63.065
## pcaPb$x[, 1]  1   57.693 63.693
summary(Pb4)
## 
## Call:
## glm(formula = Pberg_pres ~ X1_emerg + mean.mud + pcaPb$x[, 1], 
##     family = "binomial", data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -0.8403     0.3504  -2.398   0.0165 *
## X1_emerg       1.0235     0.4987   2.052   0.0401 *
## mean.mud       0.6400     0.4165   1.537   0.1244  
## pcaPb$x[, 1]  -1.1073     0.6659  -1.663   0.0963 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 64.104  on 49  degrees of freedom
## Residual deviance: 54.721  on 46  degrees of freedom
## AIC: 62.721
## 
## Number of Fisher Scoring iterations: 4

We need to extract the results from comparison between models:

# Extract AIC values
aic_values <- c(AIC(Pb), AIC(Pb1), AIC(Pb2), AIC(Pb3), AIC(Pb4))
names(aic_values) <- c("Pb", "Pb1", "Pb2", "Pb3", "Pb4")

# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min

# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))

# Create a summary table
model_comparison <- data.frame(
  Model = names(aic_values),
  AIC = aic_values,
  Delta_AIC = delta_aic,
  Akaike_Weight = aic_weights
)
print(model_comparison)
##     Model      AIC Delta_AIC Akaike_Weight
## Pb     Pb 69.35803  6.637270    0.01860592
## Pb1   Pb1 67.44012  4.719367    0.04854208
## Pb2   Pb2 65.60165  2.880890    0.12171344
## Pb3   Pb3 63.81621  1.095456    0.29719446
## Pb4   Pb4 62.72076  0.000000    0.51394410
summ(Pb4)
Observations 50
Dependent variable Pberg_pres
Type Generalized linear model
Family binomial
Link logit
χ²(3) 9.38
Pseudo-R² (Cragg-Uhler) 0.24
Pseudo-R² (McFadden) 0.15
AIC 62.72
BIC 70.37
Est. S.E. z val. p
(Intercept) -0.84 0.35 -2.40 0.02
X1_emerg 1.02 0.50 2.05 0.04
mean.mud 0.64 0.42 1.54 0.12
pcaPb$x[, 1] -1.11 0.67 -1.66 0.10
Standard errors: MLE
plot_summs(Pb4, inner_ci_level = .9, colors =  "black", plot.distributions = T)

sumbb<-summary(Pb4)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Pelophylax bergeri"

coefPb<-dfcoef

coef_final<- full_join(coef_final, coefPb)
## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`

Rana temporaria

length(data$points[data$Rtemp_pres==1])
## [1] 28

PCA

pcaRt<-prcomp(data[,c(11:15)])
pcaRt
## Standard deviations (1, .., p=5):
## [1] 0.5465325 0.4866070 0.4235086 0.2588685 0.2377199
## 
## Rotation (n x k) = (5 x 5):
##                      PC1         PC2         PC3        PC4         PC5
## Bbuf_pres    0.625089286  0.69848866 -0.02523248  0.3427587 -0.05706788
## Harbo_pres  -0.003738853  0.26551360  0.11132516 -0.6440599 -0.70871860
## Lvulga_pres  0.345240883 -0.19486754 -0.87020121 -0.2881804  0.05037168
## Pridi_pres   0.355978808  0.03901502  0.36656799 -0.5993301  0.61497015
## Pberg_pres   0.602774607 -0.63412979  0.30878437  0.1595581 -0.33724709
str(pcaRt)
## List of 5
##  $ sdev    : num [1:5] 0.547 0.487 0.424 0.259 0.238
##  $ rotation: num [1:5, 1:5] 0.62509 -0.00374 0.34524 0.35598 0.60277 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:5] 0.38 0.08 0.24 0.12 0.34
##   ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:50, 1:5] 0.391 0.035 0.0573 0.0536 -0.5678 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
PCA(data[,c(11:15)])

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
pcaRt$x[,1]
##  [1]  0.39099996  0.03502115  0.05733583  0.05359698 -0.56775346 -0.56775346
##  [7] -0.56775346 -0.56775346  0.38026203 -0.22251257  0.03502115 -0.57149231
## [13]  0.05733583  1.00535132  1.01608925 -0.56775346  0.03502115 -0.56775346
## [19]  0.40257671  0.03502115  0.05733583 -0.56775346  0.38026203  1.01608925
## [25]  0.05733583 -0.56775346 -0.22251257 -0.22251257  0.05733583 -0.56775346
## [31]  0.66011044  0.03502115 -0.56775346  0.40957578  1.36133013 -0.56775346
## [37] -0.56775346 -0.56775346  0.38026203 -0.56775346  0.03502115  1.01608925
## [43] -0.56775346  0.05359698  1.00535132 -0.56775346  0.40257671  0.05733583
## [49] -0.56775346  0.40257671
Rtsp<- glm(Rtemp_pres~ pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rtsp, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2033251 0.1676742 0.6665984 0.983952 0.01463262 0.03435206 0.2463931 0.3776317 0.03524028 0.08412745 0.2870765 0.5998516 0.4891765 0.07965571 0.5448528 0.2563778 0.3082929 0.07674179 0.9848365 0.1750712 ...
summary(Rtsp)
## 
## Call:
## glm(formula = Rtemp_pres ~ pcaRt$x[, 1], family = "binomial", 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.2482     0.2916   0.851    0.395
## pcaRt$x[, 1]  -0.8049     0.5492  -1.466    0.143
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 68.593  on 49  degrees of freedom
## Residual deviance: 66.341  on 48  degrees of freedom
## AIC: 70.341
## 
## Number of Fisher Scoring iterations: 4
Rt<-glm(Rtemp_pres~X1_submerg + X1_emerg +  X1_emerg_Land +
          mean.mud + I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rt, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.332263 0.2432032 0.9342931 0.992455 0.02245004 0.06821337 0.4206712 0.3812086 0.05518246 0.1453942 0.4090205 0.5073173 0.8272215 0.1083123 0.6544879 0.4293556 0.5339506 0.06591923 0.9941302 0.180026 ...
testZeroInflation(simulateResiduals(Rt))

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.98956, p-value = 1
## alternative hypothesis: two.sided
drop1(Rt)#submerged 
## Single term deletions
## 
## Model:
## Rtemp_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + 
##     I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[, 1]
##               Df Deviance    AIC
## <none>             52.815 70.815
## X1_submerg     1   52.819 68.819
## X1_emerg       1   52.868 68.868
## X1_emerg_Land  1   53.663 69.663
## mean.mud       1   52.892 68.892
## I(mean.mud^2)  1   55.137 71.137
## forest.dist    1   55.497 71.497
## MeanTemp       1   53.199 69.199
## pcaRt$x[, 1]   1   56.994 72.994
Rt1<-glm(Rtemp_pres~X1_emerg +  X1_emerg_Land +
          mean.mud + I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rt1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3339161 0.2432032 0.9342931 0.9923353 0.02265049 0.06772263 0.4206712 0.3776317 0.05545564 0.1453942 0.4090205 0.5073173 0.8272215 0.1073409 0.6544879 0.4293556 0.5339506 0.06591923 0.9941302 0.180026 ...
drop1(Rt1)# emerg
## Single term deletions
## 
## Model:
## Rtemp_pres ~ X1_emerg + X1_emerg_Land + mean.mud + I(mean.mud^2) + 
##     forest.dist + MeanTemp + pcaRt$x[, 1]
##               Df Deviance    AIC
## <none>             52.819 68.819
## X1_emerg       1   52.869 66.869
## X1_emerg_Land  1   53.663 67.663
## mean.mud       1   52.895 66.895
## I(mean.mud^2)  1   55.139 69.139
## forest.dist    1   55.568 69.568
## MeanTemp       1   53.212 67.212
## pcaRt$x[, 1]   1   57.143 71.143
Rt2<-glm(Rtemp_pres~X1_emerg_Land +
          mean.mud + I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rt2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3355691 0.2432032 0.9342931 0.9922155 0.02265049 0.06870412 0.4296856 0.3812086 0.05545564 0.147223 0.401399 0.5298257 0.8197094 0.1083123 0.6578101 0.4262667 0.5339506 0.06099988 0.9937634 0.1833293 ...
drop1(Rt2)# mean temp
## Single term deletions
## 
## Model:
## Rtemp_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + forest.dist + 
##     MeanTemp + pcaRt$x[, 1]
##               Df Deviance    AIC
## <none>             52.869 66.869
## X1_emerg_Land  1   53.678 65.678
## mean.mud       1   52.919 64.919
## I(mean.mud^2)  1   55.183 67.183
## forest.dist    1   55.964 67.964
## MeanTemp       1   53.359 65.359
## pcaRt$x[, 1]   1   57.205 69.205
Rt3<-glm(Rtemp_pres~X1_emerg_Land +
          mean.mud + I(mean.mud^2) + forest.dist + pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rt3, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2942429 0.2235656 0.9294259 0.9929341 0.02124736 0.07655602 0.3966328 0.4742061 0.04917249 0.1572818 0.391237 0.5598368 0.793417 0.1092838 0.6411988 0.4262667 0.5434854 0.06001601 0.9941302 0.1849809 ...
drop1(Rt3)# emerged Land 
## Single term deletions
## 
## Model:
## Rtemp_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + forest.dist + 
##     pcaRt$x[, 1]
##               Df Deviance    AIC
## <none>             53.359 65.359
## X1_emerg_Land  1   54.550 64.550
## mean.mud       1   53.756 63.756
## I(mean.mud^2)  1   55.441 65.441
## forest.dist    1   57.349 67.349
## pcaRt$x[, 1]   1   58.037 68.037
Rt4<-glm(Rtemp_pres~mean.mud + I(mean.mud^2) + forest.dist + pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rt4, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2231617 0.2265868 0.8393832 0.9911377 0.01603575 0.05005586 0.3876185 0.4026695 0.05272383 0.153624 0.2972385 0.5798442 0.6882474 0.1058838 0.6677769 0.4262667 0.3082929 0.06395149 0.9955977 0.2378325 ...
drop1(Rt4)# mean mud 
## Single term deletions
## 
## Model:
## Rtemp_pres ~ mean.mud + I(mean.mud^2) + forest.dist + pcaRt$x[, 
##     1]
##               Df Deviance    AIC
## <none>             54.550 64.550
## mean.mud       1   56.119 64.119
## I(mean.mud^2)  1   56.257 64.257
## forest.dist    1   60.775 68.775
## pcaRt$x[, 1]   1   59.729 67.729
Rt5<-glm(Rtemp_pres~forest.dist + pcaRt$x[,1], 
        data, family = "binomial")
simulateResiduals(Rt5, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2578758 0.2145022 0.7761099 0.9881436 0.01463262 0.0471114 0.3816089 0.3812086 0.04534796 0.1316778 0.248969 0.6973878 0.6206384 0.08985553 0.6544879 0.3428667 0.2606187 0.1023224 0.9910731 0.2692132 ...
drop1(Rt5)
## Single term deletions
## 
## Model:
## Rtemp_pres ~ forest.dist + pcaRt$x[, 1]
##              Df Deviance    AIC
## <none>            60.542 66.542
## forest.dist   1   66.341 70.341
## pcaRt$x[, 1]  1   62.849 66.849
summary(Rt5)
## 
## Call:
## glm(formula = Rtemp_pres ~ forest.dist + pcaRt$x[, 1], family = "binomial", 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    0.3043     0.3139   0.969   0.3323  
## forest.dist    0.7675     0.3426   2.240   0.0251 *
## pcaRt$x[, 1]  -0.8772     0.5937  -1.478   0.1395  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 68.593  on 49  degrees of freedom
## Residual deviance: 60.542  on 47  degrees of freedom
## AIC: 66.542
## 
## Number of Fisher Scoring iterations: 4

We need to extract the results from comparison between models:

# Extract AIC values
aic_values <- c(AIC(Rt), AIC(Rt1), AIC(Rt2), AIC(Rt3), AIC(Rt4), AIC(Rt5))
names(aic_values) <- c("Rt", "Rt1", "Rt2", "Rt3", "Rt4", "Rt5")

# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min

# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))

# Create a summary table
model_comparison <- data.frame(
  Model = names(aic_values),
  AIC = aic_values,
  Delta_AIC = delta_aic,
  Akaike_Weight = aic_weights
)
print(model_comparison)
##     Model      AIC Delta_AIC Akaike_Weight
## Rt     Rt 70.81472 6.2649552    0.01735929
## Rt1   Rt1 68.81935 4.2695814    0.04707843
## Rt2   Rt2 66.86926 2.3194892    0.12481853
## Rt3   Rt3 65.35867 0.8089063    0.26564276
## Rt4   Rt4 64.54977 0.0000000    0.39806112
## Rt5   Rt5 66.54157 1.9918035    0.14703987
summ(Rt5)
Observations 50
Dependent variable Rtemp_pres
Type Generalized linear model
Family binomial
Link logit
χ²(2) 8.05
Pseudo-R² (Cragg-Uhler) 0.20
Pseudo-R² (McFadden) 0.12
AIC 66.54
BIC 72.28
Est. S.E. z val. p
(Intercept) 0.30 0.31 0.97 0.33
forest.dist 0.77 0.34 2.24 0.03
pcaRt$x[, 1] -0.88 0.59 -1.48 0.14
Standard errors: MLE
plot_summs(Rt5, inner_ci_level = .95, colors = "black", plot.distributions = T)

sumbb<-summary(Rt5)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Rana temporaria"

coefPr<-dfcoef

coef_final<- full_join(coef_final, coefPr)
## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`

Figure 4

fig.4 <- ggplot(coef_final, 
                aes(factor(variable, 
                           levels = c("(Intercept)", "I(forest.dist^2)", "forest.dist", "I(mean.mud^2)", "mean.mud", "MeanTemp", "X1_emerg", "X1_emerg_Land", "pcaPb$x[, 1]", "pcaPr$x[, 1]", "pcaRt$x[, 1]"), 
                           labels = c("Intercept", "Distance to Forest^2", "Distance to Forest", "Mud Depth^2", "Mud Depth", "Mean T°C", "Emerged Vegetation", "Emerged Land", "PCA P. bergeri", "PCA P. ridibundus", "PCA R. temporaria")),
                    Estimate)) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="black") +
  geom_errorbar(aes(ymin = Estimate - `Std. Error`, ymax = Estimate + `Std. Error`, colour = species_name), lwd = 1, width = 0) +
  geom_point(size=3, aes(colour = species_name)) +
  scale_color_manual(values = palette2) +
  coord_flip() + 
  facet_grid(.~species_name) +
  guides(colour = "none") +
  labs(x="", y="Estimates") + 
  theme_minimal() + 
  theme(strip.text = element_text(face = "italic"))

fig.4